BUCKLING OF GRADED COATINGS 


A CONTINUUM MODEL 


Fazil Erdogan 
and 

Tz-ChengChiu 


Lehigh University 
Bethlehem, PA 18015 


October 2000 


AIR FORCE OFFICE OF SCIENTIFIC RESEARCH 
GRANT F49620-98-1-0028 


DISTRIBUTION STATEMENT A 

Approved for Public Release 
Distribution Unlimited 


20010419 101 


DTIC CONVERSATION RECORD FOR DISTRIBUTION STATEMENT REQUEST 


DTIC Personnel Making Call 

(SAOC /2/kjZ- 


Date .. Time „ , 

n oi I 


Authorizing Official 


Phone 


Agency 


/tycc*. . ******»^ */K* cLah%j£ 


Title 


Internet Document URL (if applicable) 


Distribution Statement (Please check one box) 

[^^DISTRIBUTION STATEMENT A: Approved for public release. Distribution is unlimited. 

□ DISTRIBUTION STATEMENT B: Distribution authorized to U.S. Government Agencies only. 

□ DISTRIBUTION STATEMENT C: Distribution authorized to U.S. Government Agencies and 
their contractors. 

□ DISTRIBUTION STATEMENT D: Distribution authorized to U.S. Department of Defense (DoD) 
and U.S DoD contractors only. 

□ DISTRIBUTION STATEMENT E: Distribution authorized to U.S. Department of Defense (DoD) 
components only. 

□ DISTRIBUTION STATEMENT F: Further dissemination only as directed by the controlling DoD 
office indicated below or by higher authority. 

□ DISTRIBUTION STATEMENT X: Distribution authorized to U.S. Government agencies and 
private individuals or enterprises eligible to obtain export-controlled technical data in 
accordance with DoD Directive 5230.25, Withholding of Unclassified Technical Data from 
Public Disclosure, 6 Nov 84. 


Reason for the above identified distribution statement (in accordance with DoD Directive 5230.24) 


Controlling Office 


Date of Distribution Statement Determination 


AQ Number (For DTIC-OCA Use Only) 





















BUCKLING OF GRADED COATINGS 
A CONTINUUM MODEL 


Fazil Erdogan 
and 

Tz-Cheng Chiu 


Lehigh University 
Bethlehem, PA 18015 


October 2000 


AIR FORCE OFFICE OF SCIENTIFIC RESEARCH 
GRANT F49620-98-1-0028 



Abstract 


Requirements for the protection of hot section components in many high temperature 
applications such as earth-to-orbit winged planes and advanced turbine systems have led 
to the application of thermal barrier coatings (TBCs) that utilize ceramic coatings on metal 
substrates. An alternative concept to homogeneous ceramic coatings is the functionally 
graded materials (FGM) in which the composition of the coating is intentionally graded to 
improve the bonding strength and to reduce the magnitude of the residual and thermal 
stresses. A widely observed failure mode in such layered systems is known to be interface 
cracking that leads to spallation fracture. In most cases, the final stage of the failure 
process for a thin coating appears to be due to buckling instability under thermally or 
mechanically induced compressive stress. The debonding and spallation problems are also 
observed in other layered material systems such as surface coatings in electronic devices 
and fiber reinforced composite laminates. In the case of TBCs, stress concentrations due 
to interface asperities and voids and relatively low toughness of interfaces eventually lead 
to the formation of interface cracks or generally highly weakened interfacial zones. 
Buckling instability may then occur under sufficiently high compressive stresses. 

The objective of this study is to develop a solution to the buckling instability problem 
by using continuum elasticity rather than a structural mechanics approach. The emphasis in 
the solution will be on the investigation of the effect of material inhomogeneity in graded 
coatings on the instability load, the postbuckling behavior, and fracture mechanics 
parameters such as the stress intensity factors and strain energy release rate. In plane strain 
and axisymmetric problems for an interface crack, the linear elastic small deformation 
theory gives only a trivial solution and does not produce instability regardless of the 
relative dimensions and the load amplitude. In this analysis, a nonlinear continuum theory 
is employed to examine the interface crack problem. First by using a perturbation 
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technique, the plane strain problem of a graded coating bonded to a homogeneous 
substrate containing an interface crack is reduced to an eigenvalue problem and the 
instability load is evaluated analytically. It is assumed that the applied load is a uniform in¬ 
plane compressive strain away from the crack region. The analytical solution of the 
instability problem permits the study of the effect of material inhomogeneity upon the 
inception of buckling and establishes benchmark results for the numerical solutions of 
related problems. 

To study the postbuckling behavior and to calculate the stress intensity factors and 
strain energy release rate a geometrically nonlinear finite element procedure with enriched 
crack-tip element is developed. Both plane strain and axisymmetric interface crack 
problems in TBCs with either homogeneous or graded coating are then considered by 
using the finite element procedure. It is assumed that the applied load is a uniform 
temperature drop. Comparison of the results with that obtained from the plate 
approximation shows that because of the higher constraints the plate theory predicts 
greater instability strains and lower strain energy release rates. It is also observed that 
compared with a homogeneous coating the graded coating gives lower strain energy 
release rate because of the lower thermal residual stress and higher bending stiffness. The 
failure of the coating may be examined either by comparing the calculated strain energy 
release rate with the mode dependent material toughness or by using a maximum stress- 
based rupture theory in the debonded part of the coating. 
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Chapter 1 
Introduction 


1.1 Introduction 

In order to meet the demanding performance requirements established by emerging 
technologies such as the aerospace plane, high performance engine, and microelectronics, 
it is often necessary that components perform multiple functions or exhibit characteristics 
not attainable by any single engineering material currently available. Examples include the 
blades in gas turbine engines which requires high heat and corrosion resistance as well as 
low heat conductivity and high mechanical toughness. Similarly, in machine tools, gears, 
and bearings high toughness and high wear resistance are required in the same 
components. In these cases, considerable efforts have been devoted to develop processing 
composites in which dissimilar materials are combined to take advantage of favorable 
properties of each material. The composites are usually particulate, fiber or filament 
reinforced, or layered in structure. Many of the laminated materials, thin films, and 
coatings fall into the latter category. 

One of the applications of layered composites under extensive development during the 
last few decades is using ceramic coatings to protect metallic components from high 
temperature and corrosive environments. These ceramic coatings, commonly referred to as 
thermal barrier coatings (TBCs), can lower the surface temperature of metallic 
components such as combustor heat shields, blades, and vanes in land-based and aero gas 
turbines operating with exhaust-gas and smoke-gas at temperatures around 1300°C. 
Lowering the temperature of such components not only prolongs their life but also saves a 
considerable amount of the energy required for their internal cooling. Thus, the TBC 
technology is considered by many as a viable near-term solution for the development of 
more efficient aircraft engines and stationary gas turbines. 
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A state-of-the-art TBC system used in aircraft engines typically consists of an 125-250 
fjtm thick thermally insulating ceramic layer, and a 50-125 /rm thick metallic bond coat 
between the ceramic layer and the metal component. Y 2 O 3 partially stabilized Zr0 2 (PSZ) 
is most commonly used as the ceramic layer because of its low thermal conductivity and 
relatively high coefficient of thermal expansion for a ceramic. The role of metallic bond 
coat (e.g., NiCrAlY) is to provide a good ceramic-metal bond and protect the substrate 
from high temperature oxidation since the transport of oxygen through the PSZ layers is 
relatively fast at typical turbine operating temperatures. The protection of the metallic 
substrate from oxidation is achieved by forming a thermally grown oxide (TGO) layer 
(mainly AI 2 O 3 ) at high temperature along the ceramic/bond coat interface as the oxygen 
diffusion barrier. 

A common feature of layered composites such as TBCs is that they consist of bonded 
dissimilar homogeneous materials. Because of the discontinuities in material properties 
across the interface, generally higher thermal and mechanical stresses would develop in the 
bonded structure. Consequently, the composite becomes very susceptible to cracking and 
debonding. A novel concept which is used to overcome the shortcomings of discontinuous 
material properties is to introduce an interfacial region with graded thermomechanical 
properties [1-5]. This is achieved by varying volume fractions of the constituents between 
zero and one hundred percent, thereby obtaining a continuous through-thickness material 
property and possibly microstructure variation. These inhomogeneous particulate 
composites with varying volume fractions are referred to as Functionally Graded Materials 
(FGMs). By using FGMs, it is possible to obtain both smoother stress distribution [ 6 ] and 
higher bonding strength [7], For example, in the TBC application it was shown in [ 6 ] that 
by replacing the homogeneous ceramic coating with a graded metal/ceramic layer, the 
stress singularity at the point of intersection of the free end and interface between the 
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homogeneous coating and the substrate disappears and stress distribution becomes 
considerably smoother. 

The synthesis of FGMs can be achieved through various techniques such as plasma 
spraying [8], chemical and physical vapor deposition [9], centrifugal casting [10] and 
combustion sintering [11]. By controlling both composition profile and microstructure, the 
concept of FGMs could provide great flexibility in material design and their potential for 
applications appears to be nearly unlimited. Besides the thermal barrier coatings on high 
temperature components, other current and potential applications of the concept of FGM 
includes wear-resistant coatings on load transfer components, armors or shields with 
improved impact resistance, and thermoelectric cells. 

Aside from the fabricated FGMs, there are also other materials in which physical 
properties vary continuously across the interfacial region or through the thickness as a 
natural consequence of material processing. For example, in many diffusion bonded 
materials the atomic composition of the two species vary continuously across the nominal 
interface resulting in a region with graded properties. In soil mechanics, because of the 
change in overburden pressure, the stiffness of the medium in macro scale varies in depth 
direction. There are also geological structures such as shale-sandstone interfacial zones 
which are naturally inhomogeneous. Other examples for natural FGMs are such biological 
materials as bones and sea shells. 

One major technical issue involved in developing the next generation of protective 
coatings is to assess the structural reliability of these coatings. Therefore, it is important to 
identify the relevant failure mechanisms, develop the necessary models, and perform the 
appropriate failure analysis. A major mode of failure in thermal barrier and a variety of 
other coatings is known to be cracking that leads to spallation. From the mechanics 
perspective the modes of failure may depend on the stress state in the coating. Usually 
residual tensile stresses tend to induce coating fracture normal to the interfaces. For 
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example, surface cracking usually occurs during the heating process and/or when the 
coating is under in-plane tensile loading. These cracks tend to branch into T-shaped and 
propagate along or parallel to the interface. Consequently, the coating may either peel off 
directly or spall off due to the connection of these T-shaped cracks. On the other hand, in 
the presence of interface cracks the residual compressive stresses induced by, for example, 
thermal expansion mismatch upon cooling during thermal cycling and/or bond coat 
oxidation provide the driving force for coating to blister and eventually spall off. These 
alternative failure processes exhibit very different dependences on the mechanical 
properties of the coatings and interfaces, and should be considered separately. While the 
fracture associated with tensile stresses is equally important, the focus of this study will be 
on the spallation associated with residual compressive stresses. 

1.2 Spallation of Coatings 

In the presence of an interface crack and in-plane compressive stress the final stage of 
the spallation appears to be through the buckling mechanism. From Figure 1.1 it may be 
seen that for flat plates the interface crack does not change the stress state without 
buckling because the prior normal and shear stresses at the interface are zero and the in¬ 
plane compression is unchanged and hence, there is no driving force. However, the driving 
force for crack growth at the interface arises when elastic buckling of the coating initiates. 
Buckling of the coating occurs when the in-plane compression in the coating exceeds a 
critical value. Propagation of the interface crack would take place once the crack driving 
force exceeds the fracture toughness of the material element within the vicinity of the 
crack tip. The spallation may also be caused by a brittle rupture mechanism due to 
bending. 
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(b) Crack with buckling 


Figure 1.1 A schematic showing (a) an interface crack without coating buckling 
and (b) the stress concentration that accompanies coating buckling. 

A prerequisite for the buckling mechanism associated with compressive stresses is the 
preexisting separation along or parallel to the interface. For example, an interface crack in 
a Rene-41/PSZ TBC system [12] requires a length of about 20 times coating thickness for 
buckling to occur after 800°C cooling. It had been suggested that the requisite separation 
may occur by localized void coalescence [13] at elevated temperatures. However, the 
source of inhomogeneous void formation and coalescence is still ambiguous. One likely 
source of void formation, analogous to void nucleation at grain boundaries during creep, 
is the presence of out-of-plane tensile stresses at the interface [14]. For example, the 
relatively large local tensions that develop in the presence of boundary waviness, 
especially in the presence of interface sliding, may exceed the critical void nucleation stress 
and permit the formation of stable voids. Another likely mechanism for crack initiation 

also involves interface asperities. Usually the bond coat-top coat interface is uneven with 
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local peaks and valleys rather than flat, especially in the case of thermal spraying process 
for the purpose of better bonding between bond coat and top coat. Oxidation preference is 
another source for interface asperities. For example, it is observed in the Pt-Al bond 
coat/electron beam physical vapor deposition (EBPVD) PSZ top coat TBC system that a 
thermally grown oxide (TGO) layer has wedge-like intrusions into the bond coat grain 
boundaries (bond coat thickness ~ 30 fim, grain size ~ 75 to 100 /zm) [15]. Regardless of 
the oxidation mechanism stresses would develop during the oxidation of curved surfaces 
and the sign (tension or compression) of the stresses is strongly mechanism and curvature 
sensitive [16]. In most cases the stress magnitude increases rapidly as the radius of 
curvature of the surface decreases. Furthermore, it has been shown that large tensile stress 
concentrations exist at the peak or the valley of asperities [13]. Hence, when the asperities 
satisfy certain size requirements, localized brittle cracking may occur. These cracks or 
voids may further connect with each other to form a large separation along or near the 
interface. 

One major player in the spallation of TBCs is the thermally grown oxide. Earlier 
experimental results [17-19] showed that under thermal cycling TBC spalled off during the 
cool-down period. Furthermore, the spallation occurred either at the bond coat-TGO 
interface, in TGO layer, at the TBC-TGO interface, or in the TBC layer and near the 
TBC-TGO interface. Intuitively, the TGO layer may be attributed as the weak cleavage 
plane for crack growth. Still, to understand the influence of TGO on the spallation of 
TBCs a more detailed investigation on the mechanisms of oxide growth and its material 
properties is required. 

Usually the oxide layer of thickness less than 1 /zm grows from bond coat during the 
TBC deposition process, and it continues to grow during the heat treatment ( ~ 1000°C) 
and during service stage (900-1300°C for turbine vane and blade). Many high temperature 
isothermal and thermal cycling oxidation experiments have been performed for various 
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TBC systems [17-22]. By using X-ray diffraction (XRD) these experiments had identified 
that the major oxide constituent is rhombohedral CK-AI 2 O 3 . Some other minor oxide 
phases, including Cr 203 , CoO, NiO, and various kinds of Co(Cr,Al )204 or Ni(Cr,Al )204 
spinels, were also detected for the case of MCrAlY (M being a metal) bond coat. It was 
observed from the fractograph of the surface of oxide scale after the spallation of top coat 
occurred that columnar AI 2 O 3 grains prevail as a subscale in contact with the bond coat 
and that NiO and Ni(Cr,Al )204 grows on the surface of the rather uniform AI 2 O 3 [18- 
20,23). Furthermore, such forms of oxide were observed on the spalling plane in the 
experiments. Consequently, it may be conjectured that these minor oxide phases would, 
first weaken the bonding strength of the interface and, second by work as interface 
asperities that enhance the crack initiation. 

From the location and microstructure of oxide layer one may presume that the 
mechanism of oxide growth is via the inward diffusion of oxygen through the ceramic top 
coat and the outward diffusion of cation (mainly Al 3+ ) toward the bond coat surface. 
Furthermore, investigations of self diffusion in AI 2 O 3 [24] have shown that oxygen 
transport along grain boundaries is considerably faster than within the grains. For 
aluminum in AI 2 O 3 , diffusion through both paths are roughly equally fast and about the 
same as that for oxygen transport at grain boundaries. As a result the new aluminum oxide 
tends to form at the columnar grain boundaries of the existing AI 2 O 3 . 
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Figure 1.2 Residual compression induced by internal oxide formation. 
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Oxide scales on, for example [25], FeCrAl alloys often show pronounced lateral 
growth which is again explained by the fresh oxide formation at grain boundaries. 
Consequently, residual stresses are likely to occur for such oxide structure (Figure 1.2). A 
consequence of the residual stress would be local sliding and debonding when the stress 
exceeds a critical value [13]. Note that the compressive residual stress is not likely to be 
uniform, but rather vary with position through the TGO layer such that the spatial 
variation is dictated by the relative rates of oxygen and aluminum diffusion to internal 
grain boundaries. Thickening of the oxide coating would, in this instance, occur by 
diffusive creep, involving a compression induced flux of aluminum and oxygen ions from 
the grain boundary to the bond coat-TGO interface and/or TGO-top coat interface. 

Because of the extended exposure to elevated temperature the TBC system may 
undergo significant stress relaxation. It is well known that most polycrystalline brittle 
materials are not amenable to plastic deformation by dislocation motion. Instead, for 
AI 2 O 3 in excess 1000°C, the primary mode of deformation involves diffusion. For 
columnar type microstructure subjected to biaxial compression, all of the grain boundaries 
are under essentially the same stress and atom transport occur from these boundaries onto 
the bond coat-TGO and TGO-top coat interfaces. The deformation rate is limited by the 
motion of the slowest diffusing species, along either the grain boundaries or the interfaces. 
The macroscopic deformation due to this diffusion mechanism may be characterized by 
Coble creep [25] which is essentially governed by linear viscoelastic constitutive laws. 
Note that usually the oxidation process occurs at temperatures that permit limited stress 
relaxation. At higher temperature the relaxation becomes more prevalent compared with 
the oxidation process. The stress effects are thus more likely to obtain at lower oxidation 
temperatures. 

Aside from the intrinsic residual stresses resulting from the oxidation process, the 
differences between the thermal expansion coefficients of the metal substrate, bond coat, 
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oxide scale, and top coat also lead to the development of residual stress during cooling 
down from the elevated temperature in the thermal-cycling process. Usually the thermal 
expansion coefficients are much smaller for ceramic TBC and TGO compared with 
metallic bond coat and substrate. It is also known that except from the compressive 
intrinsic stress in the TGO layer arisen from oxidation process, the TBC system is under 
nearly stress-free condition at elevated temperature. Consequently, during the cool down 
process in-plane residual compressive stresses are developed in the top coat and TGO 
layer. The residual compressive stresses may lead to interface separation and eventually 
the spallation of coating from the substrate. 

One may observe that compressive residual stresses would still be developed in the top 
coat during the cool-down process and thus the spallation may occur even without the 
existence of TGO layer. However, due to the fact that compared to other materials used in 
the TBC system the thermal expansion coefficient is the smallest and the stiffness is the 
highest for aluminum oxide (e.g. Table 1.1 [26]), the magnitude of the residual stresses in 
the TGO layer would be very high. The higher (in magnitude) residual stresses associated 
with oxide may accelerate the crack initiation and propagation and thus lead to earlier 
spallation. 

By introducing the concept of FGM replacing the homogeneous ceramic top coat, the 
in-plane thermal residual stresses would reduce significantly. However, the advantage of 
FGM coating might be compromised once the metallic particle in the FGM oxidizes. In an 
earlier experiment [27] it was observed that, again, spallation occurred in the FGM 
coating upon the cooling down stage after extended exposure to elevated temperature. 
Metallographic and XRD analysis indicated that most of the metallic phase in the spalled 
region was oxidized. One may imagine that during service an unfavorable oxide layer is 
"inserted" into the graded material transition region. A concept to overcome this is to add 
an extra oxygen barrier graded layer such that the alternative TBC system contains multi- 
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FGM layers (thermal barrier/oxygen barrier FGM layer and oxygen barrier/metallic bond 
coat FGM layer) above the structural components [27). 


Table 1.1 Material properties of a thermal barrier coating system 
at 22, 566, and 1149°C [26]. 



Young's 

Modulus 

(GPa) 

Poisson's 

Ratio 

Coefficient of 
Thermal Expansion 
( x lO-^C' 1 ) 

Substrate 

175.8 

mm 

13.91 

(Ni-based 

150.4 


15.36 

Superalloy) 

94.1 

H 

19.52 

Bond Coat 

137.9 

0.27 

15.16 

(NiCrAlY) 

121.4 

0.27 

15.37 


93.8 

0.27 

17.48 

TGO 

386 

0.257 

6 

(A1 2 0 3 ) 

349 

0.257 

8 


311 

0.257 

8.9 

TBC 

27.6 

0.25 

10.01 

(Zr02-8wt.% 

6.9 

0.25 

11.01 

Y 2 0 3 ) 

1.84 

0.25 

12.41 


1.3 Literature Survey 

Debonding or spallation induced by in-plane compressive stress is considered to be 
one of the dominant modes of failure in thin films and coatings used, for example, in 
microelectronics and optical applications as well as in wear, corrosion, and heat resistant 
coatings. In many composite laminates delaminations are also being frequently observed 
due to either manufacturing processes or low-velocity lateral impact. Many analytical and 
experimental investigations have been performed to clarify the influence of embedded 
cracks/delaminations on the compressive strength of coatings and composite laminates 
[28-38]. Surveys and reviews of such studies have been given, e.g., by Garg [39], Abrate 
[40], and Hutchinson and Suo [41]. Most of these investigations have utilized structural 
mechanics theories such as those for beams, plates, and shells, and, therefore, employed 
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approximations on the coating thickness/debonding length ratio, the boundary conditions 
enforced at the edges of debonding, and the deformed shape of the coating. Such analyses 
provide estimates of the critical debonding size as well as other practical information. 
However, solutions using a more exact theory are preferable because they provide more 
accurate results and give a basis for assessing the accuracy of those simplified models. 

Various studies have been carried out on the instability problem of embedded cracks 
within the context of continuum mechanics. By using exact equilibrium equations Keer, 
Nemat-Nasser, and Oranratnachai [42] estimated the buckling loads for a half space or 
layer containing an array of coplanar equally spaced cracks. This study was done for a 
class of hypoelastic materials for which the Jaumann rate of the Kirchhoff stress relates to 
the deformation rate tensor by an isotropic linear relation similar to Hooke's law. Wang et 
al. [43] studied the local buckling stability problem of a half space containing a through- 
the-thickness crack subjected to in-plane compression by solving the stability equation of 
linear elasticity theory and using singular integral equation technique. The plane stain 
stability problem of a layer bonded to a half space with an interface crack was also 
investigated by Wang and Takao [44]. Using a similar approach, Madenci and coworkers 
[45-50] examined the local buckling and imperfection problems of a layer or half-space 
containing a penny shaped crack subjected to mechanical or thermal loading. 

Because of the nonlinear nature of the problem, finite element methods have been used 
extensively in investigating the general problems of buckling and delamination growth. By 
using a geometrically nonlinear finite element code based on structural mechanics Nilsson 
and Giannakopoulos [51] analyzed buckling induced delamination growth of an initial 
circular delamination for composite laminates loaded in in-plane compression. Based on 
the method given by Storakers and Andersson [52], in this analysis the strain energy 
release rates were extracted from the calculated plate resultants. Nilsson and coworkers 
[53-56] further extended the same finite element procedure to account for crack-tip 
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fracture mode separation by reducing the loading at the crack tip in a nonlinear plate to 
that of a linear beam. The stress intensity factor were then identified by two-dimensional 
split beam solution by Suo and Hutchinson [57] and Suo [58]. By using a three- 
dimensional finite element formulation with modified crack closure method [59] for 
extracting fracture mechanics parameters Whitcomb [60-62] considered various relevant 
factors in the analysis of delamination buckling and growth including the influence of 
crack tip mode dependence, general contours, and contact at delamination buckling. As 
for the studies involving TBCs, Nusier and Newaz [63] investigated the thermal buckling 
problem of a stepped disk specimen with an interface crack at the bond coat-top coat 
interface by using a two-dimensional nonlinear finite element procedure with a virtual 
crack extension method [64] for fracture mechanics parameter read-out. 

Although many crack problems involving thermal or mechanical loading have been 
studied extensively for FGMs (e.g., [65-67] for problems involving interface cracking and 
[68-70] for problems with cracks perpendicular to the interface), very limited studies have 
been done for the spallation problem associated with in-plane compressive stresses. The 
only work that the author is aware of was done by Bao and Cai [71] who studied the 
delamination cracking problem in functionally graded coating/substrate system and 
examined the effect of various factors including coating gradation, location of the crack, 
and the thickness ration of the coating/substrate system by following the split beam 
approach [58], 

1.4 Overview 

The main objectives of this study of the interface crack problem in a functionally 
graded coating/homogeneous substrate material system subjected to thermally or 
mechanically induced in-plane compressive stresses may be summarized as 
• To develop a fracture model 
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• To determine the load/crack length relation at instability 

• To obtain the solution of the crack problem 

• To examine the effect of material inhomogeneity on the crack problem. 

In Chapter 2 the problem of a FGM coating bonded to a semi-infinite homogeneous 
substrate containing an interface crack subjected to fixed-grip compression parallel to the 
free surface is formulated within the context of kinematically nonlinear plane elasticity. By 
using a perturbation technique the nonlinear equilibrium equations are reduced to a pair of 
linear partial differential equations, in which the compressive strain appears as unknown 
constant coefficient, with the displacement perturbations as unknown functions. These 
equations may further be reduced to a pair of homogeneous Cauchy-type singular integral 
equations by means of applying Fourier transformation, introducing a new set of unknown 
functions (i.e., the derivative of crack opening displacements), and substituting the 
boundary and continuity conditions. The resulting eigenvalue problem is then solved 
numerically as described in Chapter 3. The compressive strain at buckling instability and 
the phase angle are obtained from the lowest eigenvalue and corresponding eigenvector of 
the homogeneous equations. By following the same procedure, the buckling instability 
problem of a FGM layer bonded to a semi-infinite homogeneous substrate with a highly 
weakened interfacial zone, which contains a series of short cracks and is modeled by a 
large crack that covers the weakened zone with crack surface bridging, is also considered 
in Chapters 2 and 3. 

Note that the perturbation technique described in Chapter 2 can only provide the 
information at the inception of elastic instability and would not be suitable to determine 
the stresses and displacements beyond the instability point. This is done by performing the 
related nonlinear postbuckling analysis in Chapter 4. In this study it is assumed that locally 
the constitutive equations are linear and the problem is solved by using a geometrically 
nonlinear finite element program. Enriched crack tip elements are used to calculate the 
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stress intensity factors and strain energy release rate directly. At each value of compressive 
load the equilibrium is determined by using the principle of virtual work. An incremental- 
iterative procedure is used for solving the related nonlinear equations. 

In Chapter 5 the results for stability and postbuckling analysis are given for some 
specific material properties. Effect of coating/substrate thickness ratio and the initial 
curvature of the specimen on the fracture mechanics parameters are studied by using finite 
element method. A more realistic axisymmetric top coat/TGO/bond coat/substrate 
specimen with an interface crack under temperature drop is also considered. The top coat 
can be either homogeneous ceramics or metal/ceramics FGM. These results are then 
compared with solutions based on structural mechanics theory. 

Chapter 6 gives the conclusions drawn from this study. Directions for future research 
are also discussed. 
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Chapter 2 

Buckling of FGM Coating - Analytical Procedure 


2.1 Governing Equations 

As shown in Appendix A, the equilibrium equations and boundary conditions of finite 
deformation theory of elasticity, without body forces, may be expressed as 


ijj + (&jk / Ui,k),j ~ 

(2.1) 

Tij + =: Ti , 

(2.2) 


respectively, where Ui, aij, Ti, and m are the components of the displacement, second 
Piola-Kirchhoff stress, applied traction, and unit normal vector, respectively. Note that 
(2.1) and (2.2) are referenced to a fixed Cartesian coordinate system under the Lagrangian 
description of deformation. The Green-Lagrange strains are given by 

= ~(uij + Uj^i “T Uk,i'U , k,j')' (2.3) 


Note that (2.1), (2.2), and (2.3) correspond to (A21), (A15), and (A9), respectively, 
except that the left sub- and superscripts used in Appendix A are omitted here for 
simplicity. A set of linearized governing equations for the buckling analysis may be 
developed under the adjacent equilibrium concept. Assume that a critical original 
equilibrium configuration, denoted by (0), exists (at bifurcation point). The displacements, 
stresses, tractions, and strains for an adjacent (buckled) configuration may be expressed as 


(0) , * 

Ui = uy + Ui, 

(2.4a) 

(0) , * 

<rij = °ij+°ip 

(2.4b) 

Ti=T^ 0) +T*, 

(2.4c) 


17 



(2.4d) 



where denotes the Kronecker delta. Similarly, by substituting ( 2 . 4 ) into ( 2 . 2 ) and (2.3) 
and neglecting higher order terms, we may obtain 


(S ik + ufj)o] k + (T^ul^nj = T*, 

(2.7) 

4 = \ [&+«13 kj++«?>;,<] ■ 

( 2 . 8 ) 

Assuming that the displacement derivatives in the original configuration are 

1, Equations (2.6)-(2.8) may further be reduced to 

small, i.e., 

<j,j + ( a fk u lk),j = o, 

(2.9) 

(a; j + af k ) ul k )n j = T;, 

( 2 . 10 ) 

€ ij = + U li)' 

( 2 . 11 ) 


Equations (2.9)-(2.11) may then be used for studying the buckling instability problems. 
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2.2 Formulation of the Crack Problem 


y 



Figure 2.1 Functionally graded coating bonded to a homogeneous substrate 
with an interface crack subjected to uniform compressive strain. 

The plane strain problem described in Figure 2.1 for a layered system subjected to a 
uniform compressive strain eo is considered. The medium consists of a FGM coating of 
thickness h bonded to a semi-infinite substrate containing an interface crack of length 2a. 
The compressive strain eo represents the magnitude of the external load. It is assumed that 
the substrate is homogeneous with elastic constants fi\, Ki, the coating is inhomogeneous 
with elastic parameters /i 2 ^ 2 , and (12 is approximated by 

H2(y) = /*ie™ (2.12) 

where /4 is the shear modulus, /c* = 3 — Avi for plane strain, 1 /,• being the Poisson's ratio. 
The subscripts, i = 1 and 2, denote the substrate and coating, respectively. It may be seen 
that at the original equilibrium configuration we have 

JO) _ _. 

e ixx ^0 > 

<4°> =0, i = 1,2. (2.13) 

Vfjy = 0 , 
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From the generalized Hooke's law for the plane elasticity problem given by 


/K + 1' 

j XX + 

**'N 

Co 

1 

\ 

\K — 1 - 

V AC — 1 . 

JfAtyy, 

/ 3 K ^ 

)/^xx + 1 

( ft + 1 N 

)fJ'€yy, 

U- 1 / 

k /c — 1 > 


(2.14) 


&xy fA^xy 

and using (2.13) we may obtain 

e (°) 


e iyy 


a<°> = 

IXX 


VKT+rk 0 ’ 

8/ijCo 


i = 1 , 2 . 


(2.15) 


Ki + 1 ’ 

Substituting (2.13) and (2.15) into ( 2 . 9 ), we have 


da* ixx 

d°i xy 

/ 8 flieo \ 

l 52 < n 

(2.16a) 

dx 

H-- 

dy 

V Ki + 1 J 

1 dx 2 °’ 

d °*xy 

d(J iyy 

( 8/. Li€ 0 \ 

, d 2 v* 

(2.16b) 

dx 

dy 

V Ki + 1 ) 

1 dx 2 °’ 

Hooke's law the a*j — e- relations may be expressed as 


* 

//c + 1 ' 

\ * / 3 

-K\ „ 


® XX 

U-i, 

+ (k 

_ 2 J^yy’ 


a yy ~ 

i i 
*— 1 ^ 

+(f 

+1 \ » 

_ 1 )^ e yy 

(2.17) 


a xy = fie 


xy * 


In previous studies it was shown that the influence of the variation in Poisson's ratio in the 
crack problem is rather insignificant and k may be assumed to be constant throughout the 
medium. Thus, by substituting (2.11) and (2.17) into (2.16), we find the governing 
equations for the stability problem as follows: 
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(2.18a) 


. ,d 2 u? . ,<y*u7 U“v7 , <./c)u7 

( K + ^IhF + ^ + 2 dxdv + ^ K ~ ^{~dv + 


d 2 u* _ d 2 v* 


'du* 


dy 2 dxdy 


-s(^) 

V K + 1 / 


dv[) 

dy ' dx ) 
d 2 u: 


eo 


dx 2 


= 0 , 


d 2 v* d 2 v 

(* -1)4§-+(«+ 1 ) 




d 2 v* 

+ 2i£+7<(3-«) 


dy 2 9x9?/ 

+ 7t( K +1) 


9y 


d^ 

dx 

>(—) 

\« + l/ 


€0 


ay 

dx 2 


(2.18b) 


= 0 , 


where * = 1, 2, and 71 = 0 for the homogeneous substrate, 72 = 7 for the FGM coating. 
From Equations (2.10), (2.13), and (2.15) the boundary conditions may be expressed in 
the general form as 

aUjrtj = T£i, k = 1,2, (2.19) 


where rij = n y along y = h and rij = —n y along y = 0 for the FGM coating; nj = n y 
along y = 0 for the homogeneous substrate. From Equation (2.19) the boundary and 
continuity conditions for the medium can be written explicitly as follows: 

<4^,0, h) = 0, <rl xy (x, h) = 0, - 00 < x < 00 , (2.20a,b) 

0i yy ( x , - 0) = ^ 2 yy( x i + 0), <r*i xy ( x , - 0) = o* 2xy {x, +0), - oo < x < oo, (2.20c,d) 

°iyy( x , - 0) = o, a* lxy (x, - 0) = 0, - a < x < a, (2.20e,f) 

ul(x, - 0) = u 2 (x, + 0), v\{x , - 0) = v 2 {x , + 0), |x| > a. (2.20g,h) 


Therefore, the buckling instability problem may be characterized by Equations (2.18) and 
(2.20). To solve Equation (2.18), first we define the Fourier transforms of the two 
displacement deviation components, u* and v*, respectively, as 


Ui(a,y) = / u*(x,y)e lax dx, 

J -OO 

* = 1, 2. (2.21) 

/ OO 

v*(x,y)e~ lax dx, 

-OO 
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By definition, u* and v* may be expressed in terms of Ui and V r , respectively as 

1 f°° 

<{x,y) = —J U l{a ,yy-da, 

i = 1 , 2 . ( 2 . 22 ) 

1 f 00 

<(z, 2 /) = -y V i (a,y)e 1QX da, 

For the inhomogeneous layer 2, by substituting (2.21) into (2.18), it may be shown 

that 

- .( K + 1 )- 8 (^Tl) e »] Q ^ + (' t - 1 )5 + 2ia! f < 2 - 23a > 

+ 7 (k - 1)~7^ + 7 (k - l)iaV 2 = 0 
dy 

- (k - 1) - s(^j-)e„ o?Vi + (* + l)^r + (2.23b) 

dV> 

+ 7(3 - /c)ia£/ 2 + 7(« + l)- = 0 


Solving Equation (2.23) we obtain that 

4 

Ui(a,y) = y]C t (o)e”‘», 

k=l 

0 <y<h, (2.24) 

4 

V 2 (a,y) = £>*(a)C*(a)e™ 

fe=i 


where Ci,..., C 4 are unknown functions of ct, n 4 are the roots of the characteristic 
equation 



and mi ,..., m 4 are given by 
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\ - /I _ 8 (* ~ 1 )CQ i - » 

^2 1/1 / , \n 1^1 j 


A 3 = 


(* + 1) 2 


i 860 | | 

1 _ JT+T 1 


\ — _ L 8 («-l)el , | 

4 y («+i ) 2 ’ 


and r 4 are given by 


n(«) = 

r 2 (a) = 
r 3 (a) = 


— l 


a 

V 


- i— 
a 


. a 

*v 


(2.29) 


(2.30) 


r 4 (a) = -i—. 

a 

However, the boundedness condition of ttj and v\ as y goes to — oo requires that 

A 3 = A 4 = 0. (2.31) 

By substituting (2.22), (2.24), (2.26) and (2.27) into (2.11) and (2.17) we may obtain the 
stress perturbation components for the inhomogeneous coating 2 (0 < y < h) as follows: 


ief) 


<T2yy(x,y) 


j«2 

27r 


/ oo 4 r 

E< 

■00 jt =1 L 




C k t nky t iax da, (2.32a) 


C k e nky e iQX da, (2.32b) 
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(2.32c) 



a* lxy (x,y) = J(Ai + ^)A l e XlV + 2\ 2 A^ d ax da. (2.33c) 

To solve the unknowns C 4 and A\, A 2 the homogeneous equations ( 2 . 20 a-d) may 
be used to eliminate four of the six unknowns and the mixed boundary conditions ( 2 . 20 e- 
h) would give a pair of dual integral equations to determine the remaining two. 

2.3 The Integral Equations 

We define the density functions as follows: 

— [v* 2 (x, + 0 ) - v\(z, - 0 )] = /1 (*), 

- 00 < x < 00 . (2.34) 

d 

— [U2(x, +0)-t*i(*, -0)] =/ 2 (x), 

It may be seen that (2.20g,h) would be satisfied if we require 

fi(x) = 0, \x\>a, * = 1,2, (2.35) 
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(2.36) 


[ fi(x)dx = 0, i = 1, 2. 

*/ —a 


Note that (2.36) is also known as the single-valuedness condition. Physically it means that 
the displacements are single-valued for the uncracked portion along y — 0. By taking 
Fourier transform of (2.34) we have 

ia [V 2 (a, + 0) - Vi (a, - 0)] = F x (a) (2.37a) 


ia[U 2 (a, + 0) - U x (a, - 0)] = F 2 (a ) 

where by definition and from (2.35) F\ and F 2 are given by 

F i (<*) = / h ( x)e~ iax dx, 

J —a 


(2.37b) 


2(a) = [ h(x) 

J —a 


(2.38) 


e~' ax dx. 


By substituting (2.32) and (2.33) into the homogeneous equations (2.20a-d) we obtain 

|j[ iQ (l^f) +min ‘(^)] e "*‘ Ct = 0' < 239 > 


4 

i 

k =1 


y^X n k + i am k )t nkh C k = 0, 


2iaA I+iQ [l + (^)^ 2+ |:[ io (^) +m ^(^±i) 


(2.40) 


C k = 0, (2.41) 


2 4 

(^i + + 2\ 2 A 2 — + iam^C^ = 0. (2.42) 


Substituting the appropriate expressions for U's and V’s into Equation (2.37) we have 


2 . 4 

- ( jj-J Ai - \ 2 A 2 + Yj am kC k = Fi, 


(2.43) 
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— ia(Ai + A 2 — C\ — C 2 ~ C% — C 4 ) — F 2 , 


(2.44) 


It is seen that all the unknown functions A\ (a), A 2 (a) and C\ (a),..., C±{a) may be 
determined in terms of linear combinations of Fi(a) and F 2 (a) by using (2.39)-(2.44). 
The detailed expression of A\, A 2 and C x ,..., C 4 in terms of F\ and F 2 are given in 
Appendix B. 

By substituting (2.33b,c) into (2.20e,f) the boundary conditions on the crack surface 
may be expressed as follows: 

Iim - r 1 2A i* Xiy + [l + f —) 2 } A 2 t^)iae iax da = 0 a < x < a, (2.45) 
y-*-° 27 ry_oo ( V a / J J 


^Urn ~2^J \ 2Aie + [1 + (^J J A 2 e A2 * ji ae lax da = 0 , - a < x < a, (2.45) 

^Um o ^ J°° |~(Ai + ^-) A^ 1 * + 2A 2 A 2 e A ^ e ax da = 0,-a<x<a. (2.46) 

From Appendix B the expression of A\ and A 2 in terms of F\ and F 2 are given by 

A _/_ ^22 _\ /_ di2 _\ F 2 (B27) 

k dud*22 ~~ d 2 idi2/ |Q!| \d\\d 2 2 ~ d2\d\2' 

A* = ( , , <fa , , lA + f -r . in , 1 )fe. (B28) 


'_^21_ 

. dnd 22 — d 2 \d\ 2 


)S+(; 


dn \ F 2 


dnd 22 — d 2 \d\2' ice 


) ia ’ 


where dij(a), i, j = 1, 2, are given in Appendix B. Now substituting (2.38), (B27), and 
(B28) into (2.45) and (2.46) and then changing the order of integrations, we may write the 
integral equations with the density functions fi and f 2 as unknowns as follows: 

CL 2 

— f '^Kjj(x,t)fj(t)dt = 0, i = 1, 2, —a<x<a (2.47) 

7T J -a , = i 


where the kernels Ki 3 are given by 

K n = Iim i /“£[AiiWe**' + 

3 /—027_oo|a| 


(2.48) 
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K n = lim - / [A 2 i(a)e A ^ + D 122 (a)t^ y ]t iQ ^da, 

(2.49) 

K 21 = lim - / [D 211 (a)e Al3/ + ^(aje^le^-^a, 

3/-^ ® l J _ OQ 

(2.50) 

^ 22 = yl™olJ j^| [-^221 (a)e Al?/ + D 222 (a)e X2y ]e ia ^ x ~^da, 

(2.51) 

in which D{jk(a), i, j, k = 1, 2, are given in Appendix B. It is observed that the 
integrands in (2.48)-(2.51) are continuous functions of a and vanish at a = 0. Therefore, 

any singularities the kernels K{j may have must come from the asymptotic behavior of the 
integrands as |a| approaches infinity. The details of the asymptotic expansions of the 
integrands in the kernels (2.48)-(2.51) are given in Appendix C. The leading terms of 
asymptotic expansions of Dijk for |a| goes to infinity are shown to be (Appendix C) 

j^oo 2A| 

D ni = !_^ 2 > (2.52a) 

(1 + AJ 2 ) 2 

112 2A^(1-Af)’ 

(2.52b) 

nOO 1+A| 2 

121 " 1-AJ 2 ’ 

(2.52c) 

noo 1 + K 2 

122 “ 1 - At 2 ’ 

(2.52d) 

noo 1+AI 2 

211 " 1 - Aj 2 ’ 

(2.52e) 

noo 1 + -M 2 

212 " l - a; 2 ’ 

(2.52f) 

noo _ (i+Aj 2 ) 2 

221 " 2AJ(1 - AJ 2 )’ 

(2.52g) 
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(2.52h) 


where 


n~ _2A|_ 
222 2 _ ^*2 * 



A2 

l<* 


(B6,7) 


In evaluating the kernels .KV, (Equations (2.48)-(2.51)) we separate the asymptotic leading 
terms (2.52) and evaluate them independently. For the case of K\\ we may write 

K£ = lim | f j^| [5m eAl ‘ Hs + D^e^je^-^da. (2.53) 


Separating the infinite integral (2.53) into two parts, — 00 to 0 and 0 to 00 , and making a 
change of variable for the part from — 00 to 0 by letting <p = — a, after some 
simplifications we obtain 

poo 

K£= lim / \D]° n e x ' ay + D£ 2 e x * ay ] sina(f - x)da. (2.54) 

y-'-QJo 


By using the relation 


poo 

lim / e, aX i y sma(t — x)da = lim 


t — x 


0 (KvY + (t- x ) 2 


(2.55) 


t — x 


, » = 1 , 2 , 


Equation (2.54) may be further simplified as 


r)°° -i- D 00 
^oo _ ■ u in ~T_- u n 2 


K% = 


t — x 


(2.56) 


which is a Cauchy kernel. After the separation of the asymptotic leading terms (2.52a,b) 
which gives the Cauchy kernel the remainder of the integrand in Kn is headed by a term 
no higher than — as |a| —► 00 as indicated in Appendix C. It is analytic everywhere and 
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we may eliminate the limiting process by substituting y = 0 directly under the integral. 
This then gives the Fredholm kernel kn (x, t) as follows 
1 f°° i(x 

*n CM) = - D? u + A i 2 (a) - Z^ 2 ] e ia ^W (2.57) 


Likewise, for the case of K 22 we have a Cauchy kernel (from the asymptotic leading 
terms) 


D<x 1 j-)oo 
r^oo _ -^221 ^ -^222 
^22 - , „ > 

(2.58) 

with a Fredholm kernel (for the remainder of K 22 ) expressed by 

MM) = \( pL [D 221 (a) - D? 21 + D 222 (q) - AI 2 ] e iQ(l_t W. 
z J -00 \ a \ 

(2.59) 

However, from (2.52c-f) Df 2l + Df 22 = Df u + Df n = 0, the Cauchy kernels do not 

l 

exist for the case of K i2 and K 2 \ and by substituting y = 0 directly under the integrals in 

(2.49) and (2.50) the Fredholm kernels ki 2 and k 2 \ can be expressed by 

1 

ky 2 (x,t) = — / [A 21 O 2 ) + Di 22 (a)]e la ^ x ~^da, 

““ 7-00 

(2.60) 

1 f°° 

k 2 i(x,t) = — J [An (<*) + D 2 \ 2 {a)\t' Ci ^ x ~ t ^ da. 

(2.61) 

It is observed that Aji — Aji + Aj2 — Ay i, j = 1,2, are real and are even functions 
of a (see Equation (2.52) and Appendix B). Thus (2.57), (2.59), (2.60), and (2.61) may 

be simplified as 

poo 

kn(x,t) = / [Dm(oc) - + D U 2 (a) ~ Du 2 )sinoi(t — x)da, 

J 0 

(2.62) 

poo 

k\ 2 (x,t) = / [A 21 (tt) — -Di 2 i + D\ 22 {q) — A^ 2 ] c °sa(i — x)da, 

J 0 

(2.63) 


30 



£>211 + £> 212 (a) - Df 12 ]cosa(t - x)da , (2.64) 



respectively. 

From the discussion above, it can be seen that the singular behavior of the kernels 
come solely from K\\ and K 22 and thus Equation (2.47) may be expressed by a system of 
singular integral equations (SIEs) as follows: 

2 r 

*vj_ a ^ [ (D ^ + (t^) + kij ^ fj ^ dt = °’ ( 2 . 66 ) 

i = 1,2, — a < x < a. 


With the single-valuedness condition (2.36) the SIEs will then be solved for the unknown 
density functions /1 and / 2 . It is interesting to note that Equations (2.36) and (2.62) are 
homogeneous. This implies that for nontrivial solution of f\ and / 2 , a corresponding 
critical value of eo would be obtained. Physically, this value (e 0 )cr represents the critical 
compressive strain at which buckling occurs. 


2.4 The Problem of a Weakened Interface 

Besides a single large separation at the coating/substrate interface as described in 
Section 2.2, a highly weakened interfacial zone may exist in the TBC systems. The 
weakened interfacial region consists of a series of small cracks with "weak" unbroken 
ligaments between them as shown in Figure 2.2(a). Mathematically we may consider the 
unbroken ligaments between x = — a and a acting as springs. Thus, the problem of a 
series of cracks may be modeled as a single crack of length 2a with springs connecting the 
crack surfaces as shown in Figure 2.2(b). The problem under consideration is the buckling 
instability of a graded coating bonded to a homogeneous substrate containing a weakened 
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interface under uniform compressive strain e 0 . Note that the problem of interest in this 
section is basically the same as the one described in Section 2.2 (Figure 2.1) except that 
crack surface tractions are not zero. 


y 



(a) 



(b) 

Figure 2.2 (a) Functionally graded coating bonded to a homogeneous substrate 
with a weakened interface subjected to uniform compressive strain, (b) Crack 
with surface spring model used for the weakened interface. 


Following the same derivation as in Section 2.2 the stability problem for the coating 
with a weakened interface may be described by (2.18) with the boundary conditions 
(2.20a-d,g,h). However, the crack surface stress boundary conditions (2.20e,f) have to be 
replaced by the ones considering the spring bridging effect. We assume that for very small 
value of bridging stresses the stresses are proportional to the crack opening displacements. 
Thus, the crack surface stress boundary conditions may be expressed as 

o\ yy (x, - 0) = si [v 2 (x, + 0) - v\ (x, - 0)], - a < x < a, (2.67a) 
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a\ X y{x, - 0 ) = s 2 [u 2 (x, + 0) - «i(x, - 0)], -a<x<a, 


(2.67b) 


where si and s 2 are the "spring constants". From the continuity condition (2.20g,h) we 
may write the crack opening displacements in terms of the density functions f\ and f 2 
defined by (2.34) as 


v 2 (x, + 0) - vi(x, - 0) = 

J 

/ fi (t)dt, — a < x < a, 

—a 

(2.68a) 

u 2 (x, + 0) - ui(x, - 0) = 

/ f 2 (t)dt, — a < x < a. 

J —a 

(2.68b) 


Following the same procedure described in Section 2.3 the stability problem for the 
weakened interface may then be expressed as 

f*a 2 

'It'D 1 ??, 4 - D??^ ( i- ) 4 . k;Ax. t) I fJt)dt - .<!.■ / f, (£\dt = 0. 

(2.69) 

i = 1, 2, — a < x < a. 


t r e f ( a?+a? 2 ) (nE)+ m*. oi w<«- ^ r /i m= o. 

^ J -a j—i L XS J J _ a 


with the single-valuedness conditions (2.36). The homogeneous integral equations (2.69) 
and (2.36) will be solved to determine the buckling strain (eo) CT - 
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Chapter 3 

Solution of Singular Integral Equations 
and Numerical Procedure 


3.1 Solution of the Integral Equations 

To solve the integral equations (2.36), (2.66), and (2.69), first the following 
normalized quantities are defined 


t x 

s = r=~, 
a a 


(3.1) 


fi (®) kij( r > s ) dkij(^x,t\ i , j — 1 , 2 . 


Equations (2.36), (2.66), and (2.69) are given, respectively, by 



i — 1, 2, 


(3.2) 


T £ E [<A?1 + ^)(^) + 


£(«)* = 0 , 

€ = 1 , 2 , - 1 < r < 1 , 


(3.3a, b) 


f/, E ps#+ D S)(£y+/;(«)** - a Si £n(s)ds =o, 


(3.4) 




3.1.1 The Infinite Integral 

In order to evaluate the Fredholm kernels k*j(r,s) for various values of the 
inhomogeneity parameter 7 h and normalized crack length a/h, we define 


__ N 
v M’ 


(3.5) 


giving 
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(3.6) 


poo 

k n ( r > s ) = (l7l^) (jr) B n (77)sin Q^h) (^j (s - r)r) dr], 

poo 

kn(r,s) = (| 7 |fc)(|) J o Si 2 (i?)cos[(| 7 |ft)(|)(. - r) V ]d V , (3.7) 

POO 

^ 2 i( r » s ) = (ItI h ^{] l )j o B 2 i(v) c os (\nf\h)(^j(s-r)r] dT], (3.8) 

poo 

^ 22 ( r > s ) = (M h ^{] l )j 0 B 22 (v)^ (h\h)(j^j(s-r)r] drj, (3.9) 

where from (2.62)-(2.65) 

B M = Da i (a) - Dfj-\ + D ij2 {a ) - D%, i, j = 1, 2. (3.10) 

Note that 77 = l/£ where £ is the expansion parameter used in the asymptotic analysis in 
Appendix C. 

The kernel k ^ is evaluated in such a way that the zero to infinity integral is broken 
into three integrals, i.e. 

&n( r > s ) = (l7|ft)(^) B ii{ri)--^- sin (|7|^)(^)(s -r)r^dr) 

+ Ia sin (\ n /\ h ){j l )( s - r )v]dv (3-11) 

+ Blu f Q (~) sin [(h\ h ) (|) ( s - r )v] dv |• 


The integrand of the first integral in (3.11) is bounded everywhere and integrated 
numerically; The third integral may be evaluated analytically and can be expressed as 


l (i)si n [(|7|A)(f)( S —r),]«i,= | 


(l7|fc)(^)(s-r) 

(bi*o(^)(s- r ) 


(3.12) 


Because both a and h are greater than zero. Equation (3.12) may be simplified as 
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(3.13) 




7T |s — r\ 

2 s - r ’ 


The second integral may be written as 

•Sin 


/ 


Si i(t?) 


V J 

r>oo 10 


sm 


d 7 i /l )(i )( s_r ) 77 dv 


/*oo 10 

+ / suW-E 

J A *--1 


(3.14) 


10 / D 

-Pllj 


s “[(l7l*)(r)(«-r)i?]di7 f 


10 (B • 

where ^ ( — 77 ^ ) is the 10 -term asymptotic expansions of Bn for large values of 77 and 


j=i 


V 3 


the expressions for Bn v j — 1, 2, are given in Appendix C. It may be seen that the 

10 {Bn \ 

asymptotic leading term of Bn ( 77 ) — ^ ( —4 j is of order I/ 77 11 and, therefore, we have 

j=i\ V 3 J 

J A S u sin[(| 7 |/i)(^)(s-r)77]d77~o(-^). (3.15) 


By choosing a large value for A in (3.15) such that 1/A 10 is small to any order that we 
desire, we may then neglect the second integral on the r. h. s. of (3.14) and express the 
Fredholm kernel kn as 


*Ti(r,«) = (|7l*)(f)|jf Su( 77 ) — sin[(|7|/t)(^)(s-r)77]d77 


+ 


J2 Bn if A (^) sin [(l7l^)(^)(s-0^]^ (3.16) 

j~2 A 


7T s — r _ 

+ 77-Sm 

2 s — r 


In Equation (3.16), only the first term is to be integrated numerically; The second term, 
which is expressed as a summation of integrals, may be determined in closed form. The 
closed form expressions for these integrals are shown in Appendix D; The last term is a 
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sign function and is discontinuous at s = r. From a previous study [72] it is shown that 
not separating this term from the numerically evaluated integral would lead to erroneous 
results, thereby inhibiting convergence of the final solution. 

We may use the same procedure used for to evaluate k 22 . The result may be 
expressed as 


&22( r > s ) - (M ■^ 22 ( 7 ?)- ~ sin (ItI^)(^)( s ~ r)rj^dTj 


+ 


10 noo i 

S 5 22 jj A (-)sin[(|7|/t)(^)(s-r)7?]d77 (3.17) 


+ 7 T --°221 

2 s — r 


where the first term will be integrated numerically and the second term has closed form 
expression. 

Having described the way to evaluate and k 22 , we now do the same for k{ 2 and 
k 2l . The kernel k* l2 (r, s) may be rewritten as 

k n(r,s) = B 12 (v)cos^\h)(^Y s ~ r )v]dv 

J a ^ 12 ( 17 )- ^p-|cos[(|7|fe)Q)(5-r)j;]d»; (3.18) 


+ 


+ B 


121 


J A (^) cos .d7|^)(^)(s —r)?/ drj 


The integrand of the first term in Equation (3.18) is bounded everywhere within the limits 
of integration; The second integral in (3.18) may be rewritten as 
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(3.19) 


poo ' 
JA . 


Bn(v) ~ 


B 


121 


V J 

oo 10 


COS 


- r )v] d V 


= L §(^) cos [ (l7 KD (s - r H*' 

(*00 10 / D 

*L 


COS [(|7|/l)(^)(s-r)7/]d7/, 


where 



is the 10 -term asymptotic expansions of B 12 for large values of 77 and 


the expressions for B^j, j — 1, 2, are given in Appendix C. The cosine integrals in the 
first term on the r. h. s. of (3.19) can be expressed in closed form (Appendix D). 
Following the same argument used in the case of k^, the second integral in (3.19) may be 
neglected if a sufficiently large value of A is chosen. The third integral in (3.18) has a 
closed form expression as follows: 

J A (“) cos [(N /l )(^)( s_r H d? l = “ Ci [(l7l^)(^)k-r|A], (3.20) 

where Ci, the cosine integral, is defined by 

r x cost — 1 

Ci(x) = 70 + lnrr + / - dt, x > 0, (3.21) 

Jo t 


in which 70 = 0.57721566490, is the Euler's constant. The reason that the integral (3.20) 
is separated from k{ 2 is similar to why the sign function is separated from the kernel kh. It 
is because the logarithmic term in Cij^(| 7 |/i) |s — r|Aj could not be properly 


approximated by numerical quadrature when integrated. As will be described in Section 
3.1.3, this term will be integrated in closed form. 

From the above discussions k\ 2 may be expressed as 
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*12 (*■>*) = 


5 i 2 ( t?)cos (\'Y\h)(j^)(s - r)rj dp 

10 poo 

+ Yl Bl2 jJ A \rf) cos (\'y\ h )(j l )( s - r )v dr] 


~ #12170 - #121 In 


#i 2 iln[(H^)(f)k 

(\y\h){f)(s-r)A CQst _ 1 


s - r\A\ 


Following the same argument as in k\ 2 we may write k 21 as 

*21 (»*»«) = (W*)(^){ J # 2 i(??)cos (\l\h)(^j( s - r )v}dv 


XU poo 1 

Yl B 2 l iJ A (^j) cos (W^)(^)(s-^]^ 

- #2ii7o - # 2 iiln[(| 7 |/i)(^)|s - r|Aj 


— #211 


p “(l'y|*)(f )(«—r)^4 


cost - 1 


(3.22) 


(3.23) 


It is interesting to note that, although not shown explicitly, k^ 2 (r,s) and k 21 (r, s) differ 
only by a sign, as is the case for problems that possesses symmetry with respect to the y- 


3.1.2 Nature of the Stress Singularity 

It is known in fracture mechanics that the stress field near the crack tip possesses a 
leading term that is proportional to r -p , where r is a small distance measured from the 
crack tip and p is the power of singularity which is a measure of the unboundedness of the 
stress field. It may be seen that if p is larger than one, the displacement field would be 
unbounded at the crack tip, which is unacceptable; on the other hand, if p is smaller than 
zero, r~ p maintains a positive power and the stress vanishes as r approaches zero and 
there is no stress singularity at the crack tip. Therefore, the value of p should be between 
zero and one. In this case the stress at the crack tip is unbounded but still integrable. 
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From the derivation of the singular integral equation in the previous chapter and the 
discussion earlier in this chapter, it can be seen that the terms in (3.3) and (3.4) that 
involve Fredholm kernels are bounded and analytic everywhere. The only contribution to 
the singular behavior is from the Cauchy kernel. Therefore, only the terms (in (3.3) and 
(3.4)) containing Cauchy kernels are considered in examining the singularities of the 
unknown density functions (s) and /£ (s). 

First we consider the singular integral equation (3.3a). By moving the nonsingular 
contributions to the right hand side, the equation can be expressed as 

\ f itt/s = p - 24 ) 

'Kj-i (s - r) 


where d>(r) represents all the bounded terms. To study the behavior of (s) near crack 
tips s = 1 and — 1 , we assume that 


/roo 


jijs) 

(1 + s)Pi (1 - s)P* ’ 


(3.25) 


where g\ ( s ) satisfies the Holder condition on — 1 < s < 1. Also, because f*(s) is of the 
same order as stress, the singularities pi and P 2 should satisfy the condition 
0 < Re(pi,p 2 ) < 1. Defining that 

Fi(z) = - f T^Kds, (3.26) 

7T J-i (S ~ Z) 


and substituting (3.26) into (3.25), we obtain 


7T J -1 


9i(s) 


(1 + s )? 1 (1 — S )P* (s — z) 


ds. 


(3.27) 


By following Muskhelishvili [73] the asymptotic behavior of (3.27) near the end points 1 
and — 1 can be expressed as 

( — l)e i7r?1 0 i(l)e-^ 


fiW = 


2^(1 + z) Pl sin( 7 rpi) 2 pi (1 — z)? 2 sin(7rp2) 


+ Fq(z), 


(3.28) 
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(3.29) 


where Fq(z) has singularity less than pi, po- Using the Plemelj formula, we obtain 


1 
7 T 



IM 

(s-r) 


ds= 1 -[F?(r) + Fr(r)}. 


From (3.28) and (3.29) it then follows that 

I f 1 IM gi(- l)cot(7Tj?i) _ g , i(l)cot(7rp 2 ) „ , , 
nj-i (s - r) 2^ (1 + r)?* 2(1 - r)^ 0 r 


Substituting (3.30) into (3.24) we have 

ffi( - l)cot(7rpi) _ ffi(l)cot(7rp2) 
2^(l+r)pi 2 p»( 1 -r)P2 


+ F 0 {r) = \D(r). 


(3.30) 


(3.31) 


Now multiplying both sides of (3.31) by (r + l) Pl and taking limit as r 

Pi(- ljcot^pQ _ 

2 U ’ 


1, we obtain 
(3.32) 


or 

cot(7rpi) = 0. (3.33) 

Similarly, multiplying both sides of (3.31) by (r — 1)^ and taking limit as r —> 1, we have 

cot ( 71 ^ 2 ) = 0. (3.34) 

The only acceptable roots of (3.33) and (3.34) that satisfy 0 < Refo,^) < 1 are 

Pi = P 2 = 1/2. By following the same procedure for other integral equations ((3.3b), 

(3.4)), we can conclude that the density function /£ ( s ) also has square root singularities at 
s = ± 1. In conclusion, we define 

/<*(«) = J^ s 2 ' ~ 1 < 5 < !> * = 1, 2, (3.35) 

where &(s), i = 1,2, satisfies the Holder condition on the closed interval of ( - 1 , 1 ) and 
9i( ~ 1) 7 ^ 0, gi( 1) 7 ^ 0, i = 1, 2. 
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3.1.3 Converting to Linear System 

Because the integral equations (3.2) and (3.3) can not be solved in closed form, a 
series expansion approach is used for solving the density functions. It may be seen that 
gi(s) in (3.35) is a bounded function. Hence, we may write gi in a series of orthogonal 
polynomials as follows: 

00 

»(«) = X£:/r,(s), *= 1,2, (3.36) 

3=0 


where T„(s) (n = 0,1,...) is the Chebyshev polynomial of the first kind and C*j are 
unknown constants. Substituting (3.35) and (3.36) into the single-valuedness condition 
(3.2) we have 



< = 1 , 2 . 


(3.37) 


Since the Chebyshev polynomials satisfy the following relation: 

i r 1 t j( 8 )ds = r i, 7 = o, 
TtJ-l y/l - S 2 \0, j > 0, 

we observe from Equation (3.37) that 

CS» = 0, < = 1,2. 


(3.38) 


(3.39) 


Because this study is to determine the buckling instability point, which corresponds to the 
onset of the first buckling mode, we may assume that the deformation is symmetric with 
respect to the y-axis. It is such that u* and v* (j = 1, 2) are odd and even functions of x, 
respectively. Furthermore, from (2.34) it may be seen that the density functions f\ and fo 
are odd and even functions of x, respectively. Since we have 

T»( - s) = ( - l) n T„(s), n = 0, 1, 2,..., (3.40) 


from (3.1), (3.35)-(3.36) it may be seen that 
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(3.41) 


Ci( 2 ») — C 2 ( 2 n-i) — 0> n = 1, 2, 3, 


Thus, we may rewrite (3.36) as 


^i( s ) — y>i/r 2j -i (s), 

j=i 


(3.42a) 


^2 (s) — E^T 2 ,( S ). 

j=l 


(3.42b) 


By substituting (3.35) and (3.42), the singular integral equation (3.3) may be rewritten as 


*j-ih \/l-5 2 [ 5-r 


+ A^n (r, s) ds 


,, r 1 00 

E 


C2k^2k{s) f 

— 7 ===k n (r, s)ds = 0, - 1 < r < 1. 


(3.43a) 


, # /»! oo 

-/ E 

* J-i h 


fi tTa(S , ) f £>gl + + fc 2 ‘ 2 (r, J 

v/l-s 2 e-r 22 V , y 


/ I oo 

,E 

■1 £=1 


^ , uT 2 fe-i(s) 
\/l - s 2 


(3.43b) 


/c|i(r,s)ds = 0, — 1 < r < 1. 


By using 


1 f 1 Tj(s )ds 
W-i (5 — r 


Uy-i(r), 


r-VT- 

r 


j = 0, \r\ < 1, 

j > 0, |r| < 1, 

f 

J>0, |r| > 1, 


(3.44) 


the singular integral equations (3.43) are regularized and can be expressed as 


/ilODm + ^ ) U2)^Cl*rU 2 A:- 2 (r) + — ^2 Clk f -J= kll ( r > s)ds 

k =1 n k =1 J y 1 S 2 

+ / ' /f ^2 ^1*2 ( r ’ s ) rfs = °> - 1 < r < 1. 

w *=1 -A-i V 1 - s 2 


(3.45a) 


&i 2 (r, s)ds = 0, - 1 < r < 1. 
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OO OO pi r-p / \ 

/MOgi +J?g 2 )Vc 2 t U a _ 1 (r) + ^y]Ca / , aW , fca(r,s)rf 

jfc=i ^ fe=i -'-ivl-s 


+ 


Mi 

7T 


oo /*1 

5°“/-. 


1 T 2fc -i(a) 
\/l - s 2 


(3.45b) 


k 2 i(r,s)ds = 0, — 1 < r < 1. 


where U n (s) (n = 0,1,...) is the Chebyshev polynomial of the second kind. The approach 
used for evaluating the Fredholm kernels k*j(r,s) has been discussed in Section 3.1.1. 
Note that in these kernels we separate the sign and logarithm functions which would 
introduce irregularities from the numerically evaluated integrals. The contribution of these 
functions in (3.45) can be expressed in closed form. Thus for the sign and logarithmic 
functions, we have 

/. 7^7 7^7^ = 7 Uj - (3 - 46) 


|s ' r|A ] rfs= ” rld 


1 w 


— r\as 


+ 




(3.47) 


By observing that 


L 


, -?==ln|s - r|ds = - ^(r), j > 1, 
-l V1 - s 2 7 


(3.48) 


from (3.38) and (3.47), we obtain 


Li (f) I s - r i A ] j - L ( 3 - 49 ) 


By substituting (3.16), (3.17), (3.22), (3.23), (3.46), and (3.49) into (3.45) and truncating 
the series at k = n. Equation (3.45) may be rewritten as follows: 
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/*! £>* [(Pf?, +I>S 2 )Ua- 2 (r) + (| 7 |A)(f)Sm 

1 /‘Tjt-iW- 


+ Mil 



Jfc=l 


7T. 


VT^ 


^n(r,s)ds 


k= 1 


+ Mil 



1 f 1 T 2*(«) 


fe=i 


7T. 


VT^ 


&12 (r,s)ds = 0, 


1 < r < 


MigC* (Pgj + D§ 2 )U 2i _,(r) + (| 7 |ft)(f )b 221 /I^7^|1^ 

+ A*i(WA)(^) gC 2i i J^J^=k 22 (r,s)ds 

K=1 

+ft (N a) (f) £ c « lj_ t ^=?ki < r - *)<**=°- 


1 < r < 


where 


fcn (r, s) = |£n (jj) - sin [(ItI^) (^) (s - 0*?] drj 

10 /»QO -| 

+ J A (—) sin [(i7l^) (£) (s — r )7?] G?77, 

kn{r,s) = S 12 (7?)cos[(|7|/i)^)(s -r)r^dj] 

10 /*oo -i 

£-Bi2i J A (^) cos [(l7l^)(^)(s-0 ? /]^ 


121 


/ 


■dt , 


(3.50a) 


1, 


(3.50b) 


1, 


(3.51) 


(3.52) 
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It may be seen that the unknowns in the linear system (3.50) are C\k and k = 1 
To solve this linear system, we use a collocation technique to reduce (3.50) to a system of 
linear algebraic equations in the unknowns C\k and C^k given by 




(3.55b) 


where k(j, i, j — 1 , 2, are given by (3.51)-(3.54) and the collocation points r m are chosen 
as 

(2m — 1)7T 

r m = cos-—-, m = 1 , n. (3.56) 

Note that although there is no restriction on the choice of the collocation points, usually 
the roots r m of Chebyshev polynomials (3.56) are chosen. These seems to give better 
convergence in numerical analysis. 

It may be seen that the system of linear algebraic equations (3.55) are homogeneous. 
For nontrivial solutions of C\k and C 2 k, k = the determinant of the matrix 

formed by the coefficients of C\k and C 2 k in (3.55) must be zero. From this condition, the 
critical compressive strain (eo) cr corresponding to the buckling instability point may be 
obtained. Furthermore, the crack opening shape at instability may be obtained from the 
corresponding C^'s and Cafe's. Note that from the nature of the eigenvalue problem it may 
be seen that the exact values of C\k s and C^t's are not available. Instead, the eigenvectors 
are determined within a multiplicative arbitrary constant. 

3.1.4 Crack Opening Shape 

From (2.68), (3.1), (3.35), and (3.42) we may express the crack opening displacement 
component in the y- and x-direction, respectively, as 

v 2 (x, + 0 ) - vx (x , - 0 ) = C lk [ ( 3 . 57 a) 

k =1 ^-1 V 1 - s 2 

X 

r = — 1 < r < 1, 

a 

u 2 (x, +0) - Uy(x, - 0) = af]C 2 k [ - 7 ^ = ds, (3.57b) 

k =1 7-iVl-s 2 

__ X 

a 

The integral in (3.57) can expressed in closed form as 
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- 1 < r < 1. 



Again, In the perturbation analysis the exact values of C^'s and C^'s can not be obtained, 
but only their relation as a "vector". Therefore, Equation (3.59) gives only the crack 
opening shape, rather than the real opening displacement. 


3.1.5 Stress Intensity Factors 

Note that because afj y = afj y = 0 (Equation (2.13)), we have a\ yy = a\ yy and 
ai xy = Oi xy . The mode I and II stress intensity factors for the crack geometry given in 
Figure 2.1 can then be defined and evaluated as follows: 

K\{ - a) = lim yj 27t( - o — ~x) <t\ vy {x, 0), (3.60) 

Ki(a) = lim\/27r(a: - o^Lia^O), (3.61) 

Kn( — a) = lim y/ 27r( — a — x) 0 \ xy {x, 0), (3.62) 

Kn(a) = limi/27r(a: - a)<j]; (x,0), (3.63) 

x—* 

where x lies outside the crack. Recognizing that while the stresses inside the crack give us 
the system of singular integral equations, these equations also provide the stress 
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expressions outside of the crack at y = 0, provided the density functions are known. 
Therefore, we can write 


*;„(*. 0) = f£ £ [( 05 , + og 2 )(^) + 


£(*)<**, (3.64) 


/ a 2 r X 

a E [<£& + + MM) /;(*)*. (3-65) 


or, in terms of the normalized quantities (3.1), 

,1 2 


= 7 £ £ [(^1 + r %)(^;) + *,>,<0 


0) = f£ £ [(Og, + o?p){jzi) + K,(r, >) 


fj{s)ds , (3.66) 


/*(s)ds. (3.67) 


By substituting (3.66) into (3.60), we have 


K,( - a) = lim V2 to( - 1 - r )/ ^ (r> S- + f J&Lds 

r - + - 1 t 7 T J-1 (s - r) 

+ T7 [ EM r » S )/J («)<*« } • 

^ «/-! j=l J 


(3.68) 


Recalling that the density functions are of the form of Equation (3.35) and k*j(r,s), i, 
j = 1, 2, are bounded in |s|<l , Equation (3.68) can be simplified as 

K,(-a) = Bm /»,(*>??, + OK 2 )£f ( - 1-r) £ (3.69) 

Following the discussion in Section 3.1.2 regarding the behavior of Cauchy integral and 
from Equations (3.26) and (3.28) with that pi = P 2 = 1/2, we have 


I [' fi(s) ds .. 9i(-l)e‘* /2 


if 

M-i 


(«-r) V^a + O 1 / 2 \/2 (1 - r) 1 / 2 


+ Fo(r), H>1, (3.70) 


where Fq is a bounded function. Substituting (3.70) back to (3.69), we obtain 
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Ki( -a) = MD? u + Dui)V™ 01 ( - !)■ (3-71) 

By substituting (3.42a) with n-term truncation into (3.71), the mode I stress intensity 
factor at x = — a can be expressed as 

71 

Ki( - a) = - MiD? n + D? u )^Y, C ^ C 3 - 73 ) 

j=\ 

where Dyh and Df l2 are given in (2.52). Similarly, from (3.61), (3.66), and (3.70) we 
may express mode I stress intensity factor at x — a as 

Kj(a) = - MDui + Du 2 )V^ 9i(l), (3-73) 

or 

n 

Ki(a) = - Ati (D? n + D? n )y/^Y,Cir ( 3 - 74 ) 

3= 1 

By following the same process, the mode II stress intensity factors at the crack tips are 
given as follows: 

n 

K n ( - a) = + Df 22 )^Y^C 2j , (3.75) 

j =i 

n 

Kn(a) = - ^(Pgj + D? 22 )^Yl C *r ( 3 - 76 ) 

3= 1 

It may be seen from the symmetry of the geometry and loading (Figure 2.1) that the mode 
I stress intensity factors would be the same at both crack tips; and the mode II stress 
intensity factors would be opposite in sign. This observation agrees with Equations (3.72) 
and (3.74)-(3.76). Note that since we can not obtain the exact values of Cu's and C 2 k's 
from the perturbation analysis but only their relative magnitudes, the values of stress 
intensity factors are not available. However, the phase angle which describes the mode 
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II/mode I ratio of stress intensity factors at the inception of buckling instability can be 
obtained as follows: 


ip(a) = tan 


-i (Ku{a) 



(3.77) 


where AJ and \\ are given in (B6) and (B7). 


3.2 On the Problem of the Weakened Interface 

The problem of the weakened interface may be described by (3.1), (3.2), and (3.4). It 
may be seen that compared with (3.3), only one extra term exists in the singular integral 
equation (3.4). Thus we may use the same solution technique with minor modifications to 
the system of linear equations. Again, we express the density functions and /£ by 
(3.35) and (3.42). The last integral in (3.4) may be written as: 

as i [ fi (s)ds = gSiy^Cit [ ^ for« = l; 

J ~ l k= i J-i V 1 - 5 2 

(3.77) 

as 2 [ fl(s)ds — as 2 T C 2k f — ds, for i — 2. 

J -i k=1 7-1VI-s 2 

By substituting (3.58) into (3.77), adding them into (3.50), and using the same collocation 
scheme we obtain the system of linear algebraic equations for the unknowns C\ k and C 2k 
as follows: 
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MiE^i* (D? n + Du2)U2k-2(rm ) + 


+ Mi (hW (j) it, Cit \f_ x %ffi n (r ™< s > ds 




?)*« iE c * 


T 2 A- (^m) 


(3.78a) 


+ Ml(l7^)(^) EC 2 4 / i ^== ^12( r mi s)ds 


+ asiE g u ^Vzy ^ \/! - = 0, m = l,...,n, 


^E C 2fe [o>£l + + (lT|^)(f )^22! 


'1 -r 


2 U 2 A:—1 Q*m) 


+ MM*)(f)«mECu%^ 

K=1 

+ Mi (M *0 (jr) E ClJfc “/ x ( r ™> s ) rfs 


(3.78b) 


U 2 fc_i(r m ) 
+ as22^<-2fc rr 


0, m = 1, ...,n, 


where k tJ , i, j = 1, 2, are given in (3.51)-(3.54) and the collocation points r m are chosen 
as (3.56). Again, it can be seen that (3.78) is a system of homogeneous equations. The 
eigenvalues obtained from (3.78) correspond to the buckling instability point and the 
eigenvectors give the corresponding buckling mode shape. 


3.3 On the Numerical Procedures 

In Section 3.1.3 we obtained a system of linear algebraic equation (3.55) (Equation 
(3.78) for the weakened interface problem). It may be seen that in the system of equations 
the coefficients of the unknowns involve some integrals that can not be evaluated 
analytically. Thus, we have to employ numerical integration technique to evaluate these 
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integrals. Two types of integrals, namely, the outer integrals involving the Chebyshev 
polynomials and the integrals from zero to A in the kernels k t /s (Equations (3.51)-(3.54)), 
are evaluated numerically by using Gaussian quadratures. 


3.3.1 Integrals Involving the Chebyshev Polynomials 

From (3.55) and (3.78) we may see that the integrals involving the Chebyshev 
polynomials may be expressed as 



% j 1) 2, 


n > 0. 


(3.79) 


A commonly used technique for evaluating (3.79) is the Gauss-Chebyshev integration 
formula [74]. Yet there are still certain restrictions when applying Gauss-Chebyshev 
formula to compute (3.79) [72], Therefore, the approach that we choose for evaluating 
(3.79) is to regularize the integrand and then to use the Gaussian quadratures for 
integration. This is achieved first by making a change of variables as follows: 


s = cos 6, 

T n (s) = cos (ncos -1 s) = cos ( n6), 


(3.80) 


giving 



cos(ra0) kij(cos6 - r m )dd, 


h j — 1) 2, (3.81) 


where from (3.51)-(3.54) it can be seen that k t j(r, s) = k t j(s — r). For better 
convergence of the numerical calculation, the integral from 6 = 0 to tt is broken into two 
integrals from 6 = 0 to cos -1 r m and from 6 = cos -1 r m to n, and are then computed by 
using Gaussian quadratures. 
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3.3.2 Convergence of the Approximation 

In Section 3.1.3 the singular integral equations are converted into a linear system of 
size 2 n by 2 n, where n is the number of terms in the summation that make up the density 
functions which are the unknowns in the singular integral equations, n is also the number 
of collocations points that we chose to form the system of linear equations. Theoretically, 
as the value of n increases, the accuracy of series approximation of the unknown density 
functions would be higher. Nevertheless, the numerical computation required for the 
higher order approximations increases in a quadratic fashion. Therefore, in order to obtain 
a satisfactory solution without spending too much computing time, it is important to 
examine the relation between the calculated solution and the number of terms n used in 
the approximation. To do that, we start a pilot case for the problem described in Figure 
2.1 by taking for the value of the material constants, Poisson's ratio, v = 0.3, the 
inhomogeneity parameter of the FGM coating, 7 h = — 2.3026 (p 2 (h)/pi = 0.1), and 
h/a = 0.3. The computed critical or instability strain (eo) CT as a function of n is given in 
Table 3.1. It can be seen from Table 3.1 that by using 8 approximation terms and 
collocation points, the results agree with higher order approximations up to 3 significant 
digits. Therefore, we will use n = 8 in our analysis unless otherwise stated. 


Table 3.1 Effect of the number of approximation terms on the instability point 
and phase angle for the problem described in Figure 2.1 with 7 h = — 2.3026, 
v = 0.3, and h/a = 0.3. 


Approximation terms n 

Instability load (eo) CT 

Phase angle ip (degree) 

2 

0.04027 

-37.82 

4 

0.03672 

-42.15 

8 

0.03671 

-42.36 

16 

0.03671 

-42.37 

32 

0.03671 

-42.39 
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Chapter 4 

Postbuckling Analysis 


4.1 Introduction 

Although the buckling instability problem can be solved analytically by using a 
perturbation technique and the bifurcation point may be obtained as described in the 
previous chapter, the solution for the problem of FGM coating/homogeneous substrate 
containing an interface crack under compressive loading described in Figure 2.1 cannot be 
obtained beyond the inception of buckling by using the perturbation technique. Moreover, 
the fracture mechanics parameters for the crack problem, e.g., stress intensity factors, are 
nonzero only for the postbuckling state. Therefore, a nonlinear, postbuckling analysis 
which studies the full displacement field is required to complete the fracture analysis. In 
the present chapter this is done by using a finite element method. The finite element 
procedure developed for solving the postbuckling problem is based on a general 
displacement-based finite element code FRAC2D [75] in which special enriched crack-tip 
elements are used to calculate the fracture mechanics parameters. 

Because nonlinear behavior is involved in the postbuckling of layered material 
containing an interface crack, an incremental-iterative finite element approach is used. 
Nonlinearities come into play because displacements (but not strains) and rotations are 
large. However, it is assumed that the displacements are small around the crack tip and 
that linear elastic fracture mechanics may be employed to investigate the behavior near the 
interface crack tip. The finite element analysis uses special enriched elements for the crack 
tip region as well as a geometrically nonlinear formulation for the rest of the 
coating/substrate composite to investigate the interface crack behavior simultaneously 
with the postbuckling deformation. The derivation of the finite element approach for the 
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geometrically nonlinear problem, the formulation for the special enriched elements, and 
the implementation of the material inhomogeneity are described in the following sections. 

4.2 Formulation of the Nonlinear Problem 



Figure 4.1 A general elastic medium under loading. 


Consider an arbitrary elastic medium subjected to external loading as shown in Figure 
4.1. By using the principle of virtual displacement [76], we may express the equilibrium 
equation with respect to the configuration of the medium at time t in terms of the 
following variational form 



(4.1) 


where f V is the volume of the medium at time t, t r ij is the Cartesian component of the 
Cauchy stress tensor, 6 t e t j is the virtual strain corresponding to virtual displacements 8u t 
and 6uj, and 8 t eij is given by 


8t&ij — 


1 rd(8ui) d(8uj)' 

2 l d t Xj d l Xi . ’ 


(4.2) 


t x l being the Cartesian coordinates of material point at time t, and the external virtual 
work t R is given by 
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(4.3) 


'R= [ 7f6M'V + f •ffSufd'S + 'ff’Suf, 

J*v j*s 

where t f B , *ff , and t ff are the Cartesian components of the applied body force, 
distributed surface traction, and concentrated external loading, respectively, at time t, 8uf 
and 8uf are 8ui evaluated, respectively on the surface *S and at the point where *jf is 
applied. Note that in this analysis, the motion of the medium is relative to a fixed Cartesian 
coordinate system. Also note that the postbuckling problem under consideration is static 
and consequently, time-independent. The terminology "time t" is used only to denote the 
deformed configuration which differs from the original, undeformed configuration. No 
inertia terms are included. 

The coordinates of an arbitrary point P in the body at time 0 are °xi, °x 2 , °x 3 ; at time 
t they are t x 1 , *x 2 , *x 3 , where the left superscript refer to the configuration of the 
medium, and therefore, we have 

t Xi = °Xi+ t u i , < = 1,2,3. (4.4) 

The notation scheme for coordinates and displacements is also used for other quantities 
(stress, body force, surface traction, etc.). In addition to that, a left subscript denotes the 
configuration with respect to which the quantity is measured. For example, denotes 
the body force at time t , but measured in configuration 0. The exception is that if the 
quantity under consideration in the same configuration in which it is also measured, the 
left subscript may not be used, e.g., t f B = \f B . 

An inherent shortcoming in the general application of (4.1) is that the configuration of 
the medium at time t is unknown. To overcome this difficulty we may rewrite (4.1) such 
that the quantities of interest are measured with respect to configuration 0 (the Lagrangian 
formulation) and the integral can be evaluated over the original, undeformed domain. For 
the Lagrangian formulation the general stain tensor we used is the Green-Lagrange strain 
tensor which is given as follows: 
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(4.5) 


0 e ij = 2^° Ui ’ j + + 0 Uk ’ 1 0 u k,j)‘ 

Note that while a differentiation is involved, the left subscript indicates the configuration 

d t X{ 

in which the coordinate is measured, for example, —. The corresponding stress 

Cj X j 

tensors we used for the Green-Lagrange strain tensor fa is the second Piola-Kirchhoff 
stress tensor Q<fij. As shown in Appendix A, the relation between Cauchy stress and the 
second Piola-Kirchhoff stress is given by 

f<pv\, , , 


"Tran — J O x m,i O x nj Q&ij- 


(4.6) 


Given that [77] 


&t&mn t ^i,m 


°T f) 


ip 


(4.7) 


and from (4.6) we may write 

[ <T, i 6,e, i d‘V= [ , 0 a ij 6 , 0 e i jd°V. (4.8) 

J*v J°v 

By assuming that the external virtual work l R defined by (4.3) is deformation-independent 
and can be specified at undeformed configuration, l R can be expressed as 

*R= \ t 0 f?6u i d°V+ [ IffSufdPS + ltfSu?. (4.9) 

J°v J°s 

Therefore, we may rewrite Equation (4.1) as follows: 

f 1 o<T ij S l 0 e i jdPV = t R, (4.10) 

J°v 

where *R is given in (4.9). 

Since in general the medium may undergo large displacements and large strains and 

therefore the kinematic relation (4.5) can not be linearized. Equation (4.10) can not be 

solved directly. Nevertheless, we may employ an incremental-iterative method to solve 

(4.10). That is, instead of applying full amount of the external loadings, a series of loading 
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increments are applied until the summation of the increments reaches the full amount. For 
each loading increment, an approximate solution can be obtained by referring all variables 
to a previously calculated known equilibrium configuration and linearizing the resulting 
equations for the increment. The solution is then improved by iterations. 

To derive the governing equations for the incremental-iterative approach, first we 
consider the equilibrium at time t + At. Similar to (4.10), we may write 

[ t+A o&i} f> t+A oCij d°V = t+At R, (4.11) 

J»v 

where 

t+At u i = t u < ? ) + u i , i = 1,2,3, (4.12) 

t+A Q<Jij = + 0 CTij, i,j= 1,2,3, (4.13) 

t+A o€ij = + oey, *, 3 = 1, 2, 3, (4.14) 

in which we assume that t uf \ V-**, and are known quantities obtained solving 
(4.10); and U{, 0 &ij, and o^ij are the incremental displacement, stress, and strain, 
respectively. For generality we also assume that the values of t uf \ and t e^ are not 
the "exact" solutions from (4.10) but with small errors. The "error" which represents the 
difference between the external virtual work CR) and the internal virtual work (1. h. s. of 
Equation (4.10)) evaluated with the approximate solutions is defined by 

*£(*) = t+M E ( o) = t R _ f t>)^>) d 0y (4. 15 ) 

Joy 

In Equation (4.15) t+At E^ represents the "intrinsic error" to be taken into account in the 
analysis at time t + At. Note that there is, again, no time history involved in the analysis 
and the terminologies "time t" and "time t + At" are only used to denote two successive 
incremental steps. By substituting (4.5) and (4.12) into (4.14) we may express the 
incremental strain 0 Cy as 
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0 tij — 0&ij + 0 Viji 


( 4 . 16 ) 


where 0 e jj and 0 Vij represent the linear and nonlinear part of 0 e^, respectively, and are 
given by 


0 e ij 2 (° Ui ’j ° U j>* 0 u k,i 0 Uk,j + 0 Uk,i 0 u lJ)i (4.17) 

1 

OVij 2 0^k,i (4.18) 

The virtual strain 8 <+A q e,^ can then be written as 

8 t+A o e ij — faij — 80 e t j + 8 oT)ij, i, j = 1 , 2, 3. (4.19) 

By substituting (4.13) and (4.19) into (4.11), we have the governing equation as 

J ov 0 on faij d°V + lSoVij d°V = t+At R - J la^ S 0 e i3 d°V. (4.20) 

It may be seen from (4.17) that 

/ ^l; 1 dicey/ ‘o4 > S‘ 0 eyd°V. (4.21) 

Joy Joy 

By substituting (4.15) and (4.21) into (4.20) we have 

f o a H fat j d °V + f faff SoVij d°V = ( t+At R - *R) + t+At E {0) . (4.22) 

joy Joy 

To linearize (4.22) we approximate the virtual strain increment and constitutive equations 
by the following relations 

8 o€ij = 8 o&ij: j = lj 2, 3, (4.23) 

0 on £ 0 C ijki 0 e k i + oCnjAT, i, j = 1,2, 3, (4.24) 
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giving 


[ oCijki oeki 6 0 eij d°V + [ 6 0 Vijd°V 
J»v J°v 

= { t+At R - t R) + t+At E {0) - f oCTijATSoeijdPV, 

J°v 


(4.25) 


where oCijki and oQrij are functions of °X{ (i = 1, 2, 3) for linear elastic inhomogeneous 
materials, t+At R - l R corresponds to the incremental external loading, AT is the 
temperature change in the medium from time t to t + At, and AT is assumed to be a 
known function of (°xi, °X 2 , °X 3 ). The linearized governing equation (4.25) is then solved 
by using finite element technique. Details of the finite element analysis will be discussed 
later. Assuming that the displacement, strain, and stress approximations have been 
calculated, we may then examine the difference between the external virtual work ( t+At R) 
and internal virtual work (1. h. s. of Equation (4.11)). Denoting the calculated values with 
a superscript (1) in anticipation that an iteration will in general be necessary, the "error" 
due to linearization is 

t+Atg(l) _ t+AtR 


/ 

J<V 


‘ +i „vS 


(4.26) 


The "error" described in (4.26) represents the "out-of-balance virtual work" resulting from 
linearization of (4.22). Also it may be seen that (4.26) has the same form as (4.15). In 
order to further reduce the "out-of-balance virtual work" we need to perform an iteration 
in which the following equation is solved until the difference between the external and 
internal virtual work is negligible within a certain convergence measure: 

l v oCijki oe<“» to!" 1 fv + h 6 °i n) (4.27) 

= n = 2, 3,..., 


where 
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(4.28) 


t+At j^(n) _ t-\-At 



and the incremental strain corrections 0 e^ and qt^ correspond to the displacement 
correction tq- which are given by 


. .J n ) — I 

° e ij 2 


(n) . 

0 u\J + 0 U 


(n) 


t-t-At, 

0 


(n-l) (n) 

4,i 0 ulj 


+ Ottfc,, 


(«) t+At 0 l (n-l) 
0 U k,j 


(4.29) 


on, 


(n) _ 


( n ) 


ij 2 ° U k,i ° U k,j ' 


(4.30) 


The displacements are then updated for each iteration as 


■«<„(»> = i+ai M jn-i) + i=1>2> 3. 


(4.31) 


Observing the resemblance between the governing equations for the incremental part 
(Equation (4.25)) and the iterative part (Equation (4.27)) of the analysis, we may further 
reduce the governing equations for the incremental-iterative approach to a single equation 
given by 


,+a 0 v<;- 1) d°v 


[ oCijuoetf SoeV d°V + / 

Jo\ 

( t+At R - t R) - f oCTt 3 AT8oe\] ] d 0 V 
Joy 


(4.32) 


+ t+At E (n- 1), n = l j2 , 3,..., 


where <5 nl is the Kronecker delta and 

—tJ*) t+At (0) _ t (*) <+At (0) _ t (*) ■ - 9 o 

u i - u i > o%- - . o e ti - 0 e i7 > 7 = 1, 2, 3, 


0 &ij O&iji OVij 0 Viji b 3 1) 2, 3. 


(4.33) 

(4.34) 


The finite element implementation of the governing equation (4.32) is described in the 
following section. 
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4.3 Finite Element Procedure for the Nonlinear Problem 


To derive the governing equations for the displacement-based isoparametric 
continuum finite elements, we consider only the incremental-iterative nodal-point 
displacements of as degrees of freedom. By substituting 




,( n ) 


° e b' “ Q » 


da r 


'*3 


d<£ 


(n) 


= - or]ij — a^8a^ 

° % dd m) dd n) t T ’ 


(4.35) 


into (4.32), from principle of virtual work we obtain 


/ 0 Ci 
J°v 


8oeL’ ) 9oe& , j0 


dot 


= 6, 


Til 


( t+At Rr 


~ t Rr)~ \ < 
J°V 


oCnjAT 


5oe 


Jr-dV+ / f+A o 4 n ~ 1} Jr^ dPv 

n) J*v J dal n) da { r ] 

(i) 

jo 


a 


(n) 


dai ^ 




q_ t+Af £!( n-1 ), 

r = 1, 2,..., M, 


(4.36) 


where M is the total number of degrees of freedom in the finite element analysis and 


f t f B d Ui n) jO V , f t f s du? (n) ,0n , t f c du ? [n) 

lv ofi ^F dV+ L ofi ^w dS+ofi ^ 


l Rr 


hs daX 


da? 


(4.37) 


r dit n) 

t+At Rr = / M if, B ^d°V + 
J°V dar’ 


[ t+At f s9uf^_, 0q , t+AtrcM 
Jos ofi dai n) ofi ~ 


c( n ) 


(9 Or 


(n) 


(4.38) 


<+a< £ , (”) = < 


o ( n ) 

t+Aip / *+Af («) ° e u jOta „ / n 

R ’~k ° a,i aa<"> ’ ^ ’ 


^i+A^( 0 ) 


ip / f-fAf ( 0 ) 0 U jOt/ „ pi 

* J.V ° a « 34°' ’ 


(4.39) 


It may be seen that Equation (4.36) represents a system of equations for the unknown 
nodal displacements. 


63 



For axisymmetric or plane problems we may replace the degrees of freedom a T for the 
general problem by incremental nodal point x- and y-displacement u\ and u* where k 
denotes the node number in the element. Following typical finite element formulation, we 
may write (4.36) in a matrix form for a single finite element with M nodal points as 
follows: 

( t+At 0 K ( r 1] + t+ %K^[ 1] )u {n) = 6 n \ t+At 0 F T + t+At 0 El n ~ l \ n = 1,2,(4.40) 


where U is the vector of increments in the nodal point displacements, t+A ^K^ and 
t+At 0 Kxl are the linear and nonlinear strain incremental stiffness matrices at time t + At, 
t+A qF t is the vector of increments in the equivalent loads applied at nodal points, and 
t+A lE^ is the vector of "out of balance" loads applied at nodal points. These vectors and 
matrices are given by 



juj (n) , ul (n \ u\ {n \ u 


2( n ) 
2 > 


Af ( n ) M( 

• 5 a \ ? a 2 




(4.41) 


^ At Q K^= [ t+At 0 Bp T C t+At 0 B[ n) d°V, (4.42) 

J°v 

f Bl L '+ M 0 sW B NL i°V, (4.43) 

Joy 


t+Ai 

o- 

F r = t+At R r - *R r - [ IB^ 1 0 S T d°V, 

(4.44) 



J°v 



f ,+a, R t - f < +i 'S W d°V, n*0. 


t+^ E {n) = < 

f J°v 

t T> ft r>(*) 

±t r —\ 0 ±> L 

Tt 0 S i * ) d°V, n = 0, 

(4.45) 


L J°v 


t+At R r = [ N Tt+At 0 f B d°V + 

f N sJ t+At 0 f s d°S + N cj t+At 0 f c 

(4.46) 

J Q V 


hs 


*Rr= [ 

N T if B d°V + 

f N sTt 0 f s d°S + N cTt 0 f c , 

(4.47) 

J°v 

hs 
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where t+At Q B^ and Bpji are the linear and nonlinear strain-displacement transformation 
matrices, respectively, C is the constitutive relation for the medium, and 

t+At 0 S ; are the stress matrix and vector, respectively, oSj is the thermal stress vector 
that results from the temperature change, t+At R T is the externally applied nodal point load 
at time t + At, N, N s , N c are the volume-, surface-, and point-displacement 
interpolation matrices, and t+At 0 f c , t+At 0 f s , and t+At 0 f B are the vectors of concentrated, 
surface, and body forces applied at time t + At. Note that in (4.45) the superscript (*) 
indicates that the values used for calculating the corresponding matrices are taken from 
the final results given by the iterations at a certain increment. The matrices and vectors 
used in (4.42)-(4.45) are given in Appendix E for two-dimensional axisymmetric, plane 
strain, and plane stress problems. Note that we use the approximated kinematic relation 
(4.23) and (E6) to compute the the incremental strain vector and to update the total strain 
vector as follows: 

„?<”>= (4.48) 

,+A <e= oe 1 " 1 + * +A $e <—». (4.49) 

For the two-dimensional problem under consideration the surface and line integrals 
replacing the volume and surface integrals defined in (4.42)-(4.47) are evaluated by using 
Gaussian integrations. For example, the incremental linear stiffness matrix t+At 0 K^ in 
(4.42) may be expressed as 

m m 

tW o*i” = £X>»,det (J) {j (C) (4.50) 

t=l j~ 1 

where Wi is the Gaussian weight on a local coordinates, det(J) is the determinant of the 
Jacobian transforming local to global coordinates, the subscripts i and j here designate 
(the global coordinates of) the Gaussian integration points and refer to the fact that the 
parent quantities are evaluated at these points. 
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4.4 Fracture Mechanics Parameters and Crack-tip Enriched Finite Elements 



Figure 4.2 Crack at bimaterial interface between two isotropic materials. 


For the interface crack problems in the layered material systems we assume, without 
loss in generality, that locally near the crack tip the displacements are small and that the 
linear elastic fracture mechanics may be used to describe the behavior around the crack 
tip. The asymptotic displacement expressions for an interface crack between two dissimilar 
isotropic materials (Figure 4.2) are given by [78] 

t U\j = fij l K\ + gij 

7 = 1,2, (4.51) 

tu 2j = fij *K\ + g2j *Ku, 


where the subscript j denotes the material number, fj and g t j are shown in Appendix F, 
t K\ and t Kn are the modes I and II stress intensity factors and they are defined such that 
*Ki + i*Kn = lim r 1/2_i£ (^ ee + i la Td \ (4.52) 

T — >0 

where the bimaterial constant e is given by 



(4.53) 


66 



The phase angle which may be used to describe the relative ratio of mode II/mode I 
contributions is defined as 


■0 = tari 



The strain energy release rate for the crack is given by 

G -2^r ) (w I + ^‘ K ' + ' K ^ 


(4.54) 


(4.55) 


where 

I Ei , for plane stress problems, 

E ( 4 - 56 ) 

- x — , for plane stain and axisymmetric problems. 

* 

The enriched element method, first introduced by Benzley [79] to treat crack problems 
in isotropic elastic media, includes the effects of the singularity (see (4.51) and (4.52)) in 
the element by "enriching" the conventional finite element displacement formulations with 
terms that give the proper singularity at the singular node. It may be seen that the 
displacements within the enriched elements can be expressed in terms of the nodal 

displacements and and stress intensity factors *K\ and l Ku as follows: 

M / M \ / M \ 

E N * %+' K i {hi - E n ^j + ‘ k 4 an - E - 

i \ 1 / v k =i / 

(4.57) 

M / M \ / M 

‘U2, = E N ” % + ' K I [hi - E N *fii ) + ’ k A an - E htAi 

k =1 \ fe=l / \ fc=l 



where M and iV* are the number of nodes in the element and the interpolation functions, 
respectively (see Appendix E). The unknowns in (4.57) are the 2 M nodal displacements 
and the modes I and II stress intensity factors. For example, for a 12-noded quadrilateral 
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element there are 26 unknowns. This type of element is sometimes referred to as a 
subparametric element. 



Figure 4.3 Schematic of element arrangement around a crack tip, 
E: enriched element, T: transition element. 


In order to integrate the enriched element formulation with the geometrically nonlinear 
finite elements described in previous section, an incremental-iterative form of the 
governing equations for the enriched elements are required. Due to the linear kinematic 
relationship in the Lagrangian formulation for the linear elastic theory, the governing 
equations for the enriched elements may be obtained by simply replacing the "full" 
quantities, e.g., total nodal displacements and stress intensity factors, with the incremental 
or iterative quantities in the governing equations for conventional linear finite elements. 
One important issue in integrating the enriched nonlinear finite elements is the inter¬ 
element displacement compatibility conditions between these two types of elements. This 
can be overcome by introducing "transition" elements (Figure 4.3) in which the element 
displacement field is given by 
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where the zeroing function Z varies from unity along the boundaries of the transition 
element that are in contact with crack-tip enriched elements to zero along the boundaries 
between the transition and the regular nonlinear elements. Although various forms of 
zeroing function may be used [75], in this study a linear zeroing function defined in 
Appendix F is chosen. Similarly, to overcome the incompatibility in the strain- 
displacement relations between the crack-tip and nonlinear elements we define an 
"artificial" strain-displacement relation as 

0 ^ij = 2 [ 0 u i,j “b O u j,i "b (1 §Uk,i qU/cJ • (4.59) 


It may be seen that the kinematic relation (4.59) is linear along the boundaries between 
enriched and transition elements and is "fully" nonlinear along the boundaries between 
transition and nonlinear elements. It is important, however, to note that the kinematic 
relation (4.59) is introduced only for the purpose of obtaining better computational 
convergence. 

Following the procedures for the nonlinear case described in previous section, we may 
develop the incremental governing equation for the enriched and transition elements as 
follows 


t+Ai , t+At Tr( n - 


at ry-\ u ~ 
O-^EL 


O^ENL 


|8‘ ( ” ) =«,.‘ +i *o F r +» a ‘4’' 1) , 


n = 1, 2,.. (4.60) 


where 
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c? w = «r w , if, 1 "', ifj" 1 }, 

3 = [ l *^B^C‘^lBf L d 0 V, 

Joy 

■t K^( n ) _ / dT t+A< 0(11) 75 jOt/ 

— / ■‘-'ENL O* J^ENLa V, 

Joy 


t+At rs-i* 
O-H-EL 


t+At 


t+At rp _ i+Af e> ^ r> 

0 ^ r ~ ~ -tl r 


fi 


B^„S T <t>V, 


t+At pi( n ) _ 
0 ~ 


t+At 


- / ,+a iBg T ' + "§ w <jV, „*o, 

J«v 

'Rr-t , 0 B i ;l T, (> S M d 0 V, n = 0, 

Joy 


(4.61) 

(4.62) 

(4.63) 

(4.64) 

(4.65) 


f tR r \ 

%=( •• . (4.66) 

VRr) 


In equations (4.61)-(4.66) the stress matrix and vectors t+A ^S^ n \ t+At 0 S^ n \ and 0 St are 
defined the same as those in the nonlinear formulation, t R r is given by (4.47), and l R T is 
the "singular" load vector corresponding to the unknown stress intensity factors and is a 
null vector if no enriched or transition elements are on a loaded boundary. The strain- 
displacement transformation matrices t+At 0 B^l and B E nl are given in Appendix F. 

Similar to the nonlinear formulation, the integrals in Equations (4.62)-(4.66) are 
evaluated by using Gaussian quadratures. However, a higher order integration is required 
for the enriched elements because the embedded analytical solution contains a strain 
singularity at the crack tip. 


4.5 Direction of Crack Growth 

Due to the lack of symmetry with respect to the interface in the FGM 
coating/homogeneous substrate material system, a mixed-mode condition exists regardless 
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of whether the loading is mode I or mixed mode. From the viewpoint of crack driving 
forces, it is apparent that the crack will not advance along the interface under the mixed¬ 
mode conditions, unless the plane of the crack is a weak cleavage plane. To determine the 
direction of crack growth, one may examine the stress state near the crack tip. It had been 
observed in [65] that for the problem of crack along the nominal interface of 
inhomogeneous/homogeneous medium with continuous material properties the stress state 
around the crack tip has the same asymptotic behavior as that in a homogeneous medium. 
Therefore, the asymptotic stress state in the immediate neighborhood of the crack tip may 
be written in polar coordinates as follows [80]: 

a rr {r,d) = -^=cos^[.Ki(l+sin 2 ^) +iirn(jsin0 —2tan^)] +0(r 1/2 ), 

<Jee{r, 0) = —?== cos ^ (. K\ cos 2 | sin 0 s ) + 0(r 1/2 ), (4.67) 

a T e(r,0) = —i=cos^[ifisin0-.Kii(3cos0- 1)] +0(r 1/2 ), 

2y27rr 2 

where r is the distance to the crack tip and 0 is measured counterclockwise from the 
interface. Following [81], we may estimate the direction of crack propagation by the plane 
of the maximum cleavage stress Gee around the crack tip. The crack growth (or kinking) 
angle 0 q may be obtained by solving 

— aee(r, 0 O ) = 0, a 6 e(r, 0 O ) > 0. (4.68) 

It may be seen that (4.68) is equivalent to finding the direction of the principal stress and is 
the same as solving a T e(r, 0 O ) = 0. Whichever case is used, we obtain 

cos ^ [K\ sin 0 O - #n (3 cos 0 O - 1)] = 0. (4.69) 
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0 

From (4.69) it may be seen that cos y = 0 gives 0 = ± n, which corresponds to the 

plane of crack. The second term in (4.69) provides the non-trivial solution and may be 
expressed in the following form 

ta ° 2 f “ 2^ tan f “ \ = °’ a$e( ~ r ' ^ > 0 ’ (4-70) 

4.6 On the Material Inhomogeneity of FGMs and Direction of Crack 

From the derivation for the geometrically nonlinear finite elements we may see that 
material inhomogeneity would not change the form of the system of linear equations 
(4.40). The material property functions are involved only in evaluating the Gaussian 
integrations and in computing stresses from strains at the Gaussian points. Therefore, the 
materia] inhomogeneity can be considered simply by evaluating the material properties at 
the Gaussian quadrature points from the global material property variation functions. 

A limitation of the enriched finite element procedure is that the enriched crack-tip 
elements are based on the asymptotic displacement fields for homogeneous materials. 
Thus, the size of the enriched element should be small enough such that the material 
property variation in the enriched element is negligible. 

4.7 Convergence Criterion and Summary of the Finite Element Procedure 

In performing the incremental-iterative procedure, it is clear that a convergence 
criterion is required to terminate the iteration in each incremental step. A straightforward 
criterion would be to require the out-of-balance load vector E r or E r within a preset 
tolerance of the original load increment. However, a difficulty with this criterion is that the 
displacement solution does not enter the termination criterion. In order to provide some 
indication of when both the displacements and forces are near their equilibrium values, we 
use a convergence criterion in which the increment in internal energy (i.e., the work done 
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by the out-of-balance loads through the displacement increments) is compared to the initial 
internal energy increment. Convergence is assumed to be reached when the following 
relations are satisfied [82]: 



tri( n - 1 ) 

0 -C'r 



‘ 0 F r +‘ 0 Bi°> 

t-ff , t~&(°) 

r “h o "r 


(4.71) 


where in this study the tolerance is chosen as 

e E = 1 x 10“ 9 . (4.72) 

Following the discussion in previous sections, the incremental-iterative finite element 

procedure may be summarized as follows: 

(i) The total external loading and temperature change of the medium are divided into 
certain number of increments. 

(ii) For a certain incremental step, set up the matrices and vectors K E \ K^\, F r , 
(Equations (4.42)-(4.45)) for nonlinear elements and K EE , K^ NL , 1F r , E^ 
(Equations (4.62)-(4.65)) for enriched elements. 

(iii) Construct the global stiffness matrix and forcing vector for the system of equations 
based on (4.40) and (4.60), solve for the incremental nodal displacements and stress 
intensity factors, and compute the total displacements. 

(iv) The incremental strains are obtained by using (4.48) and the total strains and stresses 
are then computed. 

(v) A new set of matrices and vectors K E \ K^\, Er l \ K EE , Kenl> 316 set U P 
with the updated quantities. 

(vi) Solving incremental and total displacements, strains, stresses, and stress intensity 
factors by following the same procedures described in (iii) and (iv). 

(vii) Repeat steps (v) and (vi) until the convergence criterion (4.71) is satisfied. 

(viii) Proceed to the next incremental step and repeat steps (ii) through (vii). 
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The finite element code based on the solution procedure described above was 
developed. The elements used in this analysis are 12-noded quadrilateral and 10-noded 
triangular cubic isoparametric elements. Details on the formulation of the elements can be 
found in [83]. To test the implementation of the code, two problems with analytical 
solutions are examined in Appendix G. Because of the lack of analytical solutions 
involving finite deformation, material inhomogeneity, and crack problems we could not 
examine the implementation of all these aspects together. However, from the excellent 
agreement between the finite element results with analytical solutions (Appendix G), it is 
safe to assume that the code implementation is working properly and that the finite 
element solutions for problems with unknown analytical solutions may be considered as 
being accurate. 
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Chapter 5 

Results and Discussion 


5.1 The Instability Analysis 

By following the numerical procedure described in Chapter 3, we can solve the 
isothermal nonlinear elasticity problem as shown in Figure 2.1. The results of interest are 
the critical or instability strain, the crack opening displacement, and the phase angle. The 
critical strain represents the bifurcation point at which buckling of the coating occurs. The 
crack opening displacement shows the shape of the crack. The phase angle shows the 
propensity of direction for the crack to grow. 

In the stability analysis, some hypothetical material combinations are used for 
comparison purpose. These combinations are: 

(1) . fi2(h)/fi\ = 20.09, 7 h = 3.0, 

(2) . p2{h)/p i = 10 ,7 h = 2.3026, 

(3) . 112 (h)/fi 1 = 0.1, 7/1 = - 2.3026, 

(4) . fi2(h) jfi\ = 0.04979, 7 h = — 3.0, 

where 7 h is the measure of material inhomogeneity 0 2 (j/) = fiitxp^y), 
7 h = \n[fi 2 (h)/fii])- Also we consider the case 7 h = 0.0001 which represents extremely 
closely to a homogeneous half space containing a crack parallel to the surface. Note that it 
is not possible to let 7 h identically equal to zero in this model due to numerical difficulties. 
It was previously shown that in crack problems for FGMs the effect of the Poisson's ratio 
on the fracture mechanics parameters is not very significant [72], Therefore, in the stability 
analysis we assume that the Poisson's ratio is the same for all materials^! =u 2 = 0.3). 

In order to validate the formulation in Chapter 2 and 3, we first examine the case of 
7 h = 0.0001 and compare the results obtained for the case of 7/1 = 0 for which an 
independent formulation is given in Appendix H. The critical strain (e 0 )cr and phase angle 
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ip for both 7 h = 0.0001 and 7/1 = 0 are given in Table 5.1. Comparison of these results 
shows that they agree with each other. Therefore, we may assume that the formulation is 
correct and the results can be considered to be accurate. 

The calculated critical value of the external load e 0 as a function of the dimensionless 
length parameter h/a for 7 h = 2.3026, 0.0001, and - 2.3026 are shown in Figure 5.1. It 
can be seen that the critical strain (eo)cr increases monotonically as h/a increases, which 
could also be seen from the classical plate theory. The figure also shows that regardless of 
the nature of coating inhomogeneity, for a given h/a the homogeneous medium always 
requires a greater instability load to initiate buckling. Intuitively one would expect that the 
critical loads for 7/2, > 0 and 7/1 < 0 would bracket (e 0 )cr for 7 h = 0. The fact that this 
argument could be very misleading may be seen from Figure 5.2 which shows the 
influence of material inhomogeneity on the buckling strain (e 0 )cr for h/a = 0.05, 0 . 1 , and 
0.15. Also shown in Figure 5.2 is the critical strain obtained from the plate theory as 
described in Appendix I. It may be seen that the critical strain given by the plate theory is 
symmetric in 7 h and becomes maximum for the homogeneous medium, 7 h = 0. The 
critical strain obtained from the continuum theory shows the same trend except that the 
maximum is slightly shifted toward negative 7 / 1 . Note that for very small values of h/a as 
expected the results obtained from plate and continuum theories are indistinguishable, 
whereas for larger values of h/a the difference between the critical strains given by the 
two theories could be very significant. This can also be seen from Figure 5.3 which shows 
the difference between the critical strains computed from continuum and plate models as a 
function of hja. The difference may be attributed to the fact that at x = =f a the plate is 
assumed to have built-in ends, whereas the continuum theory imposes no such artificial 
constraints. Furthermore, due to higher constraints at the ends, the plate model predicts 
higher values for the critical or instability load compared with those computed by 
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continuum model. This can also be observed from the results shown in Figures 5.2 and 
5.3. 

Figures 5.4-5.7 show the crack opening displacements for some fixed values of coating 
inhomogeneity 7 h and the dimensionless length parameter h/a. Note that in these figures 
the opening and shearing components (^(x, + 0) — Vy(x, — 0) and 
m 2 (x, + 0 ) — u\{x, — 0 ), respectively) are given normalized with respect to the crack 
opening at x — 0, i.e., u 2 (0, + 0) — Vy (0, — 0). As mentioned in Chapter 3 because of 
the eigenvalue-problem nature of the perturbation analysis, only the shape but not the 
exact displacement of the crack opening may be obtained. It can be seen from Figures 5.4 
and 5.5 that for larger coating thickness or shorter crack length the slope of the crack 
opening is higher. Also given in Figure 5.4 is the crack opening displacement predicted by 
plate theory which is not a function of h/a (see Appendix I). As expected, the crack 
opening displacement obtained from plate model is very close to the one for the case of 
small h/a by using continuum model. Nevertheless, as the value of h/a increases, the 
difference between the crack opening displacements around the crack tip given by 
continuum and plate models is becoming larger and larger. Again, this is due to the "fixed- 
end" boundary condition used in the plate model. Figures 5.6 and 5.7 show the influence 
of material inhomogeneity on the crack opening shapes at the inception of buckling 
instability. 

Figure 5.8 shows the phase angle ip as a function of h/a for 7 h = 2.3026, 0.0001, 
and — 2.3026. It is observed that generally the phase angles are negative, which suggests 
that the mode II stress intensity factors are negative. It can further be seen that the phase 
angle ip is more negative at the instability point for the case of smaller 7 h except for small 
h/a ( h/a < 0.1). This can also be seen from Figure 5.9 which shows the phase angle ip 
vs. coating inhomogeneity 'yh for h/a — 0.3. 
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Table 5.1 Critical strain and phase angle for a homogeneous half-space 
containing a crack parallel to the surface subjected to fixed-grip 
compression as described in Figure 2.1 with 7/1 = 0 and 7 h = 0.0001. 



Critical strain (eo) CT 

Phase angle ip (degree) 


jh = 0.0001 

7 /t = 0 

7 h = 0.0001 

0 

II 

•fL 

0.05 

0.001921 

0.001921 

-39.1 

-39.1 

0.1 

0.007103 

0.007103 

-39.8 

-39.8 

0.15 

0.01467 


-40.3 

-40.3 

0.2 

0.02381 

0.02381 

-40.6 

-40.6 

0.25 

0.03388 

6.03388 



0.3 

0.04439 

0.04439 

-40.8 

-40.8 

0.35 

0.05498 

0.05498 

-40.7 

BsU 

0.4 

0.06540 

0.06540 

-40.6 


0.45 

0.07551 

0.07551 

-40.5 

-40.5 

0.5 

0.08520 

0.08520 

-40.3 

-40.3 
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h/a 

Figure 5.3 The difference between the instability strain predicted by plate approximation 
(eo)cr-pi. and that by continuum model (e 0 ) CT <ont. as a function of h/a. 



x/a 

Figure 5.4 Normalized crack opening displacements for 7 h = — 2.3026, 
6v(x, 0 ) = v 2 (x, + 0 ) - Vi ( x, - 0 ). 
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x/a 

Figure 5.5 Normalized crack opening displacements for 7/1 = — 2.3026, 
8u{x, 0 ) — U 2 (x, + 0 ) — u\ ( x, - 0 ). 



Figure 5.6 Normalized crack opening displacements for h/a = 0.3, 
8v(x, 0 ) = vi(x, + 0 ) — vi(x, - 0 ). 
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Figure 5.7 Normalized crack opening displacements for h/a = 0.3, 
6u(x, 0) = ti 2 (x, + 0) — u\ ( x, — 0). 



Figure 5.8 Phase angle ijj for the interface crack at buckling instability 
as a function of h/a. 













5.2 The Weakened Interface Model 

In the instability analysis of plane-strain weakened interface problem described in 
Figure 2.2, two real material combinations [84] are used. They are: 

(1) . Zirconia on Rene-41 (n 2 (h)/fi i = 0.6873, 7 h — - 0.3750), 

(2) . Zirconia on Ti-6A1-4V (/z 2 (/i)M = 1-294, 7 h = 0.2577), 

where the mechanical properties of the ceramics and metal constituents are given in Table 
5.2. 


Table 5.2 Material properties used in the numerical 
calculations in Section 5.2 [84]. 


Material 

E (GPa) 

V 

Rene-41 

219.7 

0.3 

Ti-6A1-4V 

116.7 

0.3 

Zirconia 

151 

0.3 


The instability strains which are calculated by using the numerical procedure described 
in Chapter 3 are given in Figures 5.10-5.15. In these results it is assumed that the opening 
and shearing components of the bridging spring constants described in Figure 2.2(b) are 
the same, i.e.. Si = s 2 = s. Figure 5.10 shows the instability strain (e 0 )cr as a function of 
dimensionless length parameter h/a for the cases of homogeneous half space containing a 
"weakened interface" zone parallel to the surface. The instability strain (e 0 ) CT vs. h/a plots 
for the weakened interface problem in a zirconia/Ti-6Al-4V FGM coating on Ti-6A1-4V 
substrate, and a zirconia/Rene-41 FGM coating on Rene-41 substrate are given in Figures 
5.11, and 5.12, respectively. The normalized spring constant associated with the unbroken 
but weak ligaments in the weakened interface region are assumed to be as/fi 1 = 0 . 001 . 
The results for as/ni = 0 which corresponds to fully separated crack (as shown in Figure 
2.1) are also presented in Figures 5.10-5.12 for comparison. It may be seen from these 
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figures that the weakened interface model gives higher instability loads compared with the 
case of fully separated case. 

The influence of the spring constant on the instability loads for the weakened interface 
problem may be seen in Figures 5.13 and 5.14 which shows the instability strains for the 
case of 7 h — 0.2577 and 7 h = — 0.3750, respectively. In these figures one may see that 
the instability strains remain rather unchanged when as/fii < 0.001; and (eo) CT increases 
sharply as the value of as/fti exceeds 0.01. The influence of asj \i\ on instability loads 
can also be seen in Figure 5.15 which shows the instability strain as a function of coating 
inhomogeneity for several fixed values of as/fj, 1 and h/a. It may be seen that as the value 
of as/Hi increases much higher instability strains are observed for the case of negative 7 h, 
and that the symmetry of the instability strain in coating inhomogeneity gradually 
disappears as 7 h increases. The higher instability strains for the case of negative 7 h may 
be attributed to the lower nominal bending stiffness of the FGM coating compared with 
that of positive 7 h. Intuitively, one may see that the coating with softening profile 
(negative 7 h) but the same bridging spring constant as/Hi requires greater driving force 
than that for stiffening profile (positive 7 h) to induce the instability and buckling. Also 
observed in Figure 5.15 is that the increase of (€o) CT for the cases of negative 7 h is more 
significant for smaller value of h/a. 
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Figure 5.10 The instability strain (eo) C r as a function of h/a for a coating bonded 
to a homogeneous half-space containing a weakened interface region. 



h/a 

Figure 5.11 The instability strain (e 0 ) CT as a function of h/a for a Zr0 2 /Ti-6A1-4V FGM 
coating bonded to a Ti-6A1-4V substrate containing a weakened interface region. 
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h/a 

Figure 5.12 The instability strain (eo)cr as a function of h/a for a Zr02/Rene-41 FGM 
coating bonded to a Rene-41 substrate containing a weakened interface region. 



asl\\, x 

Figure 5.13 The influence of the normalized bridging constant as j\i\ on the instability 
strain (eo) CT for the weakened interface, 7 h = 0.2577, h/a = 0.3. 
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as/\x j 

Figure 5.14 The influence of the normalized bridging constant as/n\ on the instability 
strain (eo) CT for the weakened interface, 7 h = - 0.3750, h/a = 0 . 3 . 



yh 

Figure 5.15 Instability strain (eo)cr as a function of coating inhomogeneity 7 h for the 
FGM coating/homogeneous substrate containing a weakened interface under fixed-grip 
compression. 
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5.3 Postbuckling Analysis 

The postbuckling problem of the debonded coating is investigated by using the 
geometrically nonlinear finite element procedure described in Chapter 4. The top coat of 
the specimen under consideration may be either a homogeneous ceramic layer or a 
smoothly graded metal/ceramic composite. It is assumed that the specimen is under either 
plane strain or axisymmetric condition. Figure 5.16(a) shows a typical finite element mesh 
and the boundary conditions (without mechanical loading) for a two-layered 
coating/substrate system. Note that by taking advantage of the symmetry of the geometry 
and loading conditions, only the right half of the specimen is modeled (for plane strain 
problem). A vertical displacement constraint is given at the outside lower comer to 
prevent rigid body motion of the specimen. As shown in Figure 5.16(b), mesh refinement 
is done around the crack tip. The size of the enriched and transition crack-tip elements are 
chosen to be around 1 /20 of the layer thickness of the coating above the interface crack. 

5.3.1 Fixed-grip Loading 

In order to obtain the fracture mechanics parameters, e.g., strain energy release rate 
and stress intensity factors, and to compare the results of the finite element analysis with 
the analytical benchmark solutions from the stability analysis, the problem considered here 
is the fixed-grip compression problem as shown in Figure 2.1. The results are given in 
Figures 5.17-5.20. In the numerical analysis the following relative dimensions are used: 
h s /h c = 30, 2 a/h c = 40, 2 l/h c = 200, where h s , h c , 2 a and 21 are the substrate 
thickness, coating thickness, crack length, and the length of the specimen, respectively. 
The results obtained by varying h s /h c and 2 l/h c showed that the dimensions used were 
sufficiently large to properly simulate the semi-infinite medium. Note that due to the 
geometry of the medium and the nature of the loading, the "preferred" load-displacement 
path, i.e., the one corresponding to buckling and crack opening, cannot be obtained in a 
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straightforward manner from the incremental-iterative technique which starts from zero 
external load. Similar to the procedure followed for the elastica problem described in 
Appendix G, first a transverse load disturbance is introduced at (x, y) = (0,0) on the 
coating to induce initial crack opening, which is then removed at a later stage of the 
incremental-loading history to obtain an equilibrium point in the postbuckling regime for 
the original problem (without the transverse load disturbance). After that the load- 
displacement path is constructed by using the incremental-iterative method which starts 
from this computed equilibrium point. Figure 5.17 shows the crack opening displacement 
at x = 0 (denoted by 6). The analytical results obtained from the stability solution are 
indicated by small circles on the e 0 axis. The figure also shows the result for the 
homogeneous medium jh c = 0 obtained from the plate theory (Appendix I) for 
comparison. The critical or instability strains (e 0 ) cr obtained from the postbuckling analysis 
are the intersection points of the postbuckling curves with the e 0 axis. The agreement 
between the critical loads obtained from the analytical and numerical methods appears to 
be quite good. Again note that for a given crack opening the plate theory requires a 
greater compressive strain e 0 than the continuum theory, which is consistent with the 
results given in Figure 5.2. The normalized strain energy release rate and stress intensity 
factors are given in Figures 5.18 and 5.19. The normalizing quantities are given by 
Ko = E s e 1 \Jnh c and Go = (1 — u^)Kq/E s where ei = 0.002. From Figure 5.18 it may 
be seen that G = 0 for €o < (eo)cr» for a given eo > (eo)cr the material inhomogeneity has 
a significant influence on the strain energy release rate, and by and large the relative trends 
conform to intuitively expected results, namely G corresponding to yh c = 0 is bracketed 
by that obtained from 7 /i c < 0 and 7 h c > 0 and G increases with decreasing 7 h c . By 
examining the results for 7 h c = 0 in Figures 5.17 and 5.18 it may be seen that the crack 
opening displacement and strain energy release rate predicted by plate theory are 
consistently smaller than that obtained from the continuum theory. This may again be 
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attributed to a higher degree of constraint in the plate model resulting from the assumption 
of "built-in" ends. From Figure 5.19 it may be seen that generally the magnitudes of the 
mode II stress intensity factors Kn are greater than that of K\. The dominance of Kn can 
also be seen from Figure 5.20 which shows the phase angle ^ as a function of compressive 
strain for the postbuckling regime. The stronger Kn implies the tendency for the crack 
kinking into the coating and spalling the coating off. However, perhaps a more practical 
approach to studying the spallation process would be either comparing the calculated 
strain energy release rates G with the mode mixity-dependent fracture toughness G c or 
analyzing the buckled part of the coating by using a maximum stress-based rupture theory. 
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(a) 



Figure 5.16 A typical finite element mesh and boundary conditions for a 
coating/substrate system containing an interface crack, (a) FE mesh for the 
right half of the medium, (b) Enlarged crack tip region, CR1 denotes the 
crack tip. 
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0.000 0.001 0.002 0.003 0.004 

So 

Figure 5.17 Crack opening displacement 8/h c at the center of crack as a function 
of compressive strain eo, the circles on the eo axis are the instability strains (co) CT 
calculated from the analytical stability solution. 



Figure 5.18 Strain energy release rate G/Go as a function of compressive 
strain e 0 , G 0 = (1 - u^)Kq/E s , K 0 — E s eiy/irfi c , e x = 0.002. 
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Figure 5.19 Modes I and II stress intensity factors K\ and Ku as functions of 
compressive strain e 0 , K 0 = E s ei \/nh c , = 0.002. 




Figure 5.20 Phase angle as a function of compressive strain cq- 
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5.3.2 Thermal Loading 

The problem considered here is the temperature change induced coating blister in the 
metal substrate/ceramic (or FGM) coating containing an interface crack. Due to the 
mismatch in coefficients of thermal expansion between the coating and substrate, 
compressive residual stress is developed in the ceramic coating (which usually has a 
smaller thermal expansion coefficient compared with metallic substrate) during the cool¬ 
down stage of the thermal cycling process. As a result the coating buckles as the 
temperature change reaches a critical value. In this part of the study we consider a layered 
composite specimen subjected to a uniform negative temperature change. It is assumed 
that the specimen is stress-free before cooling down and that the deformation is elastic. 


5.3.2.1 The Plane Strain Problem 

As an example the plane strain problem for a nickel-based Rene-41 substrate coated 
with zirconia, containing an interface crack and subjected to a uniform temperature drop 
AT is considered. The dimensions of the specimen given in Figure 5.21 [6] are h c = 2 
mm, h s = 12.5 mm, and l = 100 mm. The relevant thermomechanical properties of the 
substrate and the coating for this specific example are given in Table 5.3. Also considered 
here is the uniform cool-down problem for the medium described in Figure 5.21 in which 
the homogeneous zirconia coating is replaced by an inhomogeneous FGM coating. A 
linear variation of the thermomechanical properties in the FGM coating is assumed, i.e.. 


E(y) = 


E c + (E, - Sc)( - -- -/*‘" ~ ), 


{ E s , 0 < y < h s , 


h s < y < h s + h c , 


(5.1) 
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[ Q c + (q s ~ a c )( — + /* c —-), h s < y < h s + h c , 

a(y) = < c (5.2) 

[ a s , 0 < y < h s . 

where (E c , q c ) and (E s , a s ) are the elastic moduli and coefficients of thermal expansion of 
the pure ceramics and metallic substrate, respectively, and are given in Table 5.3. 


Table 5.3 Material properties of the plane-strain specimen 
described in Figure 5.21 [84]. 



Material 

E (GPa) 

V 

Q (10' 6o C _1 ) 

Substrate 

Rene-41 

219.7 


16.7 

Ceramics 

Zirconia 

151 

0.3 

10.0 


Figures 5.22-5.30 show the results of the postbuckling analysis for the homogeneous 
zirconia coating. The constants used for normalization in these figures are To = 100 K, 
Kq = E s a s Toy/nh^, Go — (1 - vI)Kq/E s . It may be seen from Figures 5.22, which 
shows the normalized crack opening displacement 8/h c aX x — 0 as a function of AT /To 
for some fixed values of a/h c , that instead of having an explicit buckling point, the "crack 
opening" is a smooth function of AT /To- Figure 5.23 shows the normalized energy 
release rate G/Gq versus temperature drop for various values of a/h c . Also shown in 
Figures 5.22 and 5.23 are the results obtained by using plate theory with the assumption 
that the buckled coating has "built-in" ends. By comparing the results from continuum and 
plate theories in these figures it may be seen that, again because of the higher degree of 
constraint at the ends, the plate model underestimates the buckling deflection and strain 
energy release rate. Note that the smooth transition, rather than an explicit buckling point 
in the crack opening displacement and strain energy release rate vs. AT/To plots may be 
explained by observing that the temperature change induced residual stresses would result 
in a slightly bending of the specimen, and this bending curvature serves as the imperfection 


96 












that triggers the opening of the crack. Also, because of the bending of the specimen, the 
solution technique (adding and removing a transverse load disturbance) used for the fixed- 
grip loading problem (Section 5.3.1) is not necessary for the thermal loading problem, 
unless the substrate/coating thickness ratio is very large such that the radius of the bending 
curvature is extremely high. K\ and K n shown in Figure 5.24 are the real and imaginary 
parts of the complex stress intensity factor obtained from the continuum analysis of the 
related bimaterial interface crack problem. The phase angle as a function of AT/To is 
given in Figure 5.25. It may be seen from the sign of Kn or the phase angle that, again it 
favors the crack kinking into the coating unless, as in most cases, the interface is the weak 
fracture plane. Note that the calculated stress intensity factors K\ and Ka have the 
dimension of Pa-m 1 / 2-16 , where the bimaterial constant e is defined by (4.53). In order to 
recast these quantities in a form with conventional units (e.g., Pa-m 1 / 2 ) and to have the 
result scalable, we use the following "rotated" complex stress intensity factor [85] 

K = K 1 + iK n = KL K = (K x + i K u )L ie , (5.3) 

where L is an arbitrary length quantity which, in this study, is selected to be equal to the 
coating thickness, i.e., L = h c . The new phase angle ip as shown in Figure 5.25 can be 
related to the computed phase angle ip through the expression 

ip = ip + e In L. (5.4) 

The stresses a xx , a yy , and a xy alone y/h c = 6.25 for the case of a/h c = 20 and 
AT = 1000 °C are plotted in Figures 5.26-5.28, respectively. Note that y/h c = 6.25 
corresponds to the interface when a < x < l and to the bottom surface of the debonded 
coating for 0 < x < a. The stress a xx on the top surface of the coating (y/h c = 7.25) is 
also given in Figure 5.26. From Figures 5.26-5.28 it is observed that the thermal residual 
bending stress o xx is the dominant stress components except for around the crack tip and 
where the interface intersecting free surface. It may also be seen from Figure 5.26 that on 
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the top surface of the coating high o xx appears around x = 0. This high o xx may induce 
surface crack and fracture of the coating. The stresses a xx along x = 0 and x — a for 
a/h c = 20 are given in Figures 5.29 and 5.30, respectively. The nature of these stresses is 
mainly associated with the global thermal bending except for the crack tip where the 
singular behavior dominants. 

Figures 5.31-5.38 show the results for the FGM coating described in (5.1) and (5.2). 
For comparison the results for the homogeneous zirconia coating given in Figures 5.22 
and 5.23 are also included in Figures 5.31 and 5.32. Figure 5.33 shows the normalized 
stress intensity factors as a function of AT/To for several values of a/h c . Basically the 
results for FGM coating show the same trend as those for homogeneous coatings. It is, 
however, observed that for the same temperature drop AT both the crack opening 
displacement at the center of crack and the strain energy release rate for the FGM coating 
are considerably smaller than that for the homogeneous coating. This may be attributed to 
the decrease in thermal residual stresses because of the smooth transition of 
thermomechanical properties in FGM and the increase in the effective bending stiffness of 
the coating. As a result the FGM coating is expected to be more resistant to buckling and 
spallation resulting from cool-down. 

Figures 5.34-5.36 give the stresses a xx , a yy , and o xy , respectively, along y/h c = 6.25 
for the case of a/h c = 20 and AT = 1000 °C. The stress cr xx on the top surface of the 
coating (y/h c = 7.25) is also given in Figure 5.34. The stresses a xx along x = 0 and 
x = a for a/h c = 20 are given in Figures 5.37 and 5.38, respectively. Basically the stress 
variations for the FGM coating show the same trend as that for the homogeneous ceramic 
coating (Figures 5.26-5.30). However, one important feature of the FGM coating is that 
the positive <j xx observed in the ceramic coating top surface (Figure 5.26) does not appear 
in the FGM coating (Figure 5.34); instead, the FGM coating surface is under compression 
which may prevent the initiation of surface cracks. 
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Figure 5.22 Crack opening displacement 6/h c at the center of the crack as a function 
of temperature drop A T/To for homogeneous zirconia coating, To = 100 K. 










Figure 5.23 Strain energy release rate for homogeneous zirconia coating under uniform 
temperature drop, T 0 = 100K, K 0 = E s a s To\/^h c , G 0 = (1 - v£)Kq/E s . 



Figure 5.24 Normalized stress intensity factors for homogeneous zirconia coating as a 
function of temperature drop. To = 100 K, K 0 = E s a s To^/n\, G 0 = (1 - u^)Kq/E s , 
Ki + iK li = (K l + iKT l )K e . 
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Figure 5.25 Phase angle ip as a function of temperature drop for homogeneous 
zirconia coating, ip = ip + e In h c , To = 100 K. 



Figure 5.26 Variation of normalized a xx (x, yo) along (i) the crack surface and 
the interface (yo/h c = 6.25) and (ii) the top surface (yo/h c — 7.25) of the 
homogeneous TBC described in Figure 5.21 under 1000 °C temperature drop, 
ci/he — 20, Tq = 100 K. 
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x/a 

Figure 5.27 Variation of normalized <J yy {x, yo) along the interface of the 
homogeneous TBC described in Figure 5.21 under 1000 °C temperature 
drop, y 0 /hc = 6.25, a/h c = 20, T 0 = 100 K. 



Figure 5.28 Variation of normalized a xy ( x, y 0 ) along the interface of the 
homogeneous TBC described in Figure 5.21 under 1000 °C temperature 
drop, y 0 /h c = 6.25, a/h c = 20, T 0 = 100 K. 
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y/h c 

Figure 5.29 Variation of normalized a xx along x = 0 for the homogeneous TBC 
described in Figure 5.21 under 1000 °C temperature drop, a/h c = 20, To = 100K. 



y/h c 

Figure 5.30 Variation of normalized o xx along x/a = 1 for the homogeneous TBC 
described in Figure 5.21 under 1000 °C temperature drop, a/h c = 20, To = 100 K. 
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Figure 5.31 Crack opening displacement at the center of the crack vs. temperature 
drop in FGM and homogeneous coatings. To = 100K. 



0 8 16 


AT / T 0 


Figure 5.32 Strain energy release rate for homogeneous and FGM coatings under 
uniform temperature drop, T 0 = 100 K, K 0 = E s a s T 0 y/nh' c , G 0 = (1 - v*)K$/E s . 
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Figure 5.33 Stress intensity factors for FGM coating as a function of temperature 
drop, To = 100K, K 0 = E s a s T 0 y/rtT c , G 0 = (1 - u 2 s )Kl/E s . 



Figure 5.34 Variation of normalized a xx (x, yo) along (i) the crack surface and 
the interface ( y 0 /h c = 6.25) and (ii) the top surface ( yo/h c = 7.25) of the 
FGM coating described in Figure 5.21 under 1000 °C temperature drop, 
a/h c = 20, To = 100 K. 
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x/a 

Figure 5.35 Variation of normalized a yy (x, yo) along the interface of the FGM 
coating described in Figure 5.21 under 1000 °C temperature drop, yo/h c = 6.25, 
a/hc = 20, T 0 = 100 K. 



Figure 5.36 Variation of normalized a xy (x, y 0 ) along the interface of the FGM 
coating described in Figure 5.21, yo/h c = 6.25, a/h c = 20, T 0 = 100 K. 
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y/h c 

Figure 5.37 Variation of normalized cr xx along x = 0 for the FGM coating described 
in Figure 5.21 under 1000 °C temperature drop, a/h c = 20, To = 100 K. 



Figure 5.38 Variation of normalized o xx along x/a = 1 for the FGM coating 
described in Figure 5.21 under 1000 °C temperature drop, a/h c = 20, To = 100 K. 
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5.3.2.2 The Axisymmetric Problem 

The first axisymmetric example considered here is described in Figure 5.39 [86]. The 
disk-shaped specimen with dimensions h c = 0.556 mm, h b = 0.184 mm, h s = 9.26 mm, 
and R = 40 mm is subjected to negative temperature changes. The relevant 
thermomechanical properties of the specific specimen are given in Table 5.4. Figures 5.40 
and 5.41 show crack opening displacement at r = 0 and strain energy release rate, 
respectively, obtained from both continuum and plate models for several values of afh c . 
The normalization factors used in the axisymmetric problems are given as follows: 

T 0 = 100K, K 0 = EsasToy/nlh, Go = (1 - v s )Kq/E s . (5.5) 

It may be seen from these figures that basically the continuum and plate models agree with 

each other very well except in the case of a/h c = 40. For the case of a/h c = 40 it is not 

only that significant discrepancy is observed between results from continuum and plate 

models, but also that the results from continuum theory for a/h c = 40 do not follow the 

trend observed in other cases (a/h c = 20 and 30). The rationale for the discrepancy 

between the case of a/h c = 40 and others can be seen by examining Figures 5.42 and 

5.43. Figure 5.42 shows, for the case a/h c = 40, the evolution of the deformed shapes of 

the specimen as temperature changes. It may be seen from this figure that around the 

inception of crack opening (AT ~ 200°C) the deflection of the separated part of the 

coating exhibits the shape of first buckling mode. As temperature drops more, the shape of 

the deflection gradually changes from first buckling mode to third buckling mode. Figure 

5.43 shows the deformed shapes of the specimen under 1000°C temperature drop for 

a/h c = 30 and 40. It is observed that unlike the case of a/h c = 40, the deflections for the 

case of a/h c = 30 remains the shape of first buckling mode. It is important to note that 

because of the inability of the plate model described in Appendix I to dictate higher 

buckling modes, the results based on the first buckling mode (as shown in Figure 5.40 and 

5.41 for a/h c = 40) cannot correctly predict the postbuckling behavior. Figure 5.44 
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shows the normalized stress intensity factors as functions of temperature drop. Again, 
relaxation in stress intensity factors due to higher buckling mode for the case of 
a/h c = 40 is observed. 

By comparing the strain energy release rates G obtained from plate and continuum 
models for the cases of a/h c = 30 in Figure 5.41, one may see that the strain energy 
release rate predicted by plate model is higher than that given by the continuum model. 
This observation does not agree with the one made in the plane strain case and seems to 
contradict the "higher constraint in plate model" explanation. The reason for the difference 
between results for plane and axisymmetric geometry is that in the postbuckling analysis 
the deformed shape is asymptotically obtained for the initial postbuckling regime. As a 
consequence, the plate model would overpredict the crack opening displacement and 
strain energy release rate (Appendix I). Apparently in the case of a/h c = 30 the effect of 
the higher constraining "fixed end" assumption in plate model (e.g., observed in the case 
of a/h c = 20) is out weighted by the approximation errors of the asymptotic analysis and 
thus the plate model overpredicts the results. 


Table 5.4 Material properties of the disk-shaped specimen 
described in Figure 5.39 [86]. 



Material 

E (GPa) 

V 

a (l 0 -6° c - 1 ) 

Substrate 

Ni-based Superalloy 

175.8 

0.25 

13.91 

Bond Coat 

NiCrAlZr 

137.9 

0.27 

15.16 

Top Coat 

Zr02-8wt%Y2C>3 

27.6 

0.25 

10.01 


To examining the influence of substrate thickness on the crack opening and fracture 

mechanics parameters, we fix the crack radius/coating thickness ratio as a/h c = 20 and 

vary the substrate thickness from the specimen's original h s /h c = 16.8 to h s /h c = 8.4 

and 4.2. The results are given in Figures 5.45-5.47. Basically, it may be seen from these 

figures that the crack opening displacement at the center of the crack, the strain energy 
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release rate, and the magnitude of stress intensity factors decrease slightly as h s /h c 
decreases. Nevertheless, the results obtained from plate model varies significantly as h s /h c 
changes and the difference between plate and continuum models increases as the substrate 
thickness decreases. Therefore, it may be concluded that plate model would give better 
approximation for greater value of substrate thickness. 


no 



z 



Figure 5.39 Specimen configuration for the first axisymmetric thermal buckling problem. 
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Figure 5.40 Crack opening displacement at the center of the crack for the 
homogeneous PSZ coating under uniform temperature drop AT, To = 100 K. 



Figure 5.41 Strain energy release rate as a function of temperature drop AT 
for the homogeneous PSZ coating. To = 100 K. 
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(c) AT/To = 4 



(d) AT/To = 5 

Figure 5.42 Evolution of the displacement of homogeneous ceramic coating under 
uniform temperature change, (a) AT /To = 2, (b) AT /To = 3, (c) AT /To = 4, (d) 
AT/To = 5, To = 100 K, a/h c = 40, Scale factor = 40. 
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(a) a/h c = 30 



(b) a/h c = 40 


Figure 5.43 Deformed shape of the homogeneous ceramic coating under 
uniform temperature change, (a) a/h c = 30, (b) a/h c = 40, A T/T 0 = 10. 








Figure 5.44 Normalized stress intensity factors for homogeneous zirconia coating under 
uniform temperature drop, To = 100 K, K 0 = T s a s To\AA:> G$ = (1 — u s )Kq/E s , 
k 1 + iKn = (K l + iK n )K. 



0 5 10 

A T/T 0 

Figure 5.45 Crack opening displacement at the center of the crack for the homogeneous 
PSZ coating under uniform temperature drop, To = 100 K. 
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A T/T 0 

Figure 5.46 Strain energy release rate as a function of temperature 
drop for the homogeneous PSZ coating. To = 100 K. 



Figure 5.47 Normalized stress intensity factors for homogeneous zirconia coating 
under uniform temperature drop, To = 100 K, Kq = E s a s To\fidi c , 

Go = (1 - Vs)Kl/E s , Ki + iKo = (. K , + iK n )K. 
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The second axisymmetric specimen we examined is shown in Figure 5.48 [63]. The 
dimensions of the specimen are given as follows [63]: h c = 127 pm, h t = 10 pm, 
h], = 48.2 pm, h s = 3.175 mm, and 2 R = 25.4 mm. The specimen under consideration 
originally consists of only three layers: superalloy substrate, metallic bond coat, and top 
coat that is either a homogeneous ceramic coating or a metal/ceramics FGM coating. A 
thin layer of thermally grown oxide (TGO) which consists of mainly aluminum oxide is 
grown at the bond coat-top coat interface after an extended exposure to high temperature 
oxidative environment. The material property gradation of the FGM is assumed such that 
the FGM has the same material properties as the bond coat at the bond coat-top coat 
interface, the same properties as the pure ceramics at the TBC surface, and the property 
variations in between are given as 


E(z) = E c + (E b -E c )[ 


h c -f - ht + h\) h s — z\pi 


h c 


)"■ 


v(z) = V c + (l/b-I/ c )( 


h c + hi + h\y + h s — z \ & 


r 


a(z) = a c + (a b - a c )( 


h c + hi + /ib + “ % \ ** 


he 


y, 


(5.6) 


ht + h\, + h s < z < h c + h t + h\y + h s , 


where (£/ b , a b ) and (E c , u c , a c ) are the thermomechanical properties of the metallic 
bond coat and ceramics, respectively, pi , P 2 , and p$ are the measures of inhomogeneity of 
the FGM, and (p 1; p 2 , pz) > 0. The power-law variations for the thermomechanical 
properties are defined such that the FGM is metal-rich for 0 < (pi, P 2 , Pz) < 1; a linear 
variation from 100% metal to 100% ceramics when Pi = pz = Pz = 1» and ceramic-rich 
for the case of 1 < (pi, Pz, Pz) < oo. Two limiting cases of the material gradation are that 
the coating consists of 100% metallic bond coat material if pi, pz, and pz equal to zero; 
and that as pi, P 2 , and pz approach infinity, the coating becomes pure ceramics. As an 
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example, the profile of elastic modulus described in (5.6) is plotted for various values p x in 
Figure 5.49. The relevant thermomechanical properties for the specimen described in 
Figure 5.48 are given in Table 5.5. 


Table 5.5* Material properties of the disk-shaped specimen described 
in Figure 5.48. 



Material 


V 

a (10' 6o C 1 ) 

Substrate 

Ni-based Superalloy 

■gH 

0.26 

15.38 

Bond Coat 

NiCrAlZr 

122.9 

0.27 

15.49 

TGO 

A1 2 0 3 

320 


8.8 

Ceramics 

Zr0 2 -8 wt % Y 2 0 3 (PSZ) 

11.54 

0.25 

10.91 


The results for the specimen described in Figure 5.48 under uniform temperature drop 
are given in Figures 5.50-5.77. The normalizing quantities used here are defined in (5.5), 
in which the quantities with a subscript s correspond to the properties of the substrate 
given in Table 5.5. The normalization of the complex stress intensity factors associated 
with crack between dissimilar material interface is done by following the method described 
in the plane strain cases (Equation (5.3)) and L = h c is used as the length parameter for 
the "rotation". It is assumed that the TBC system is stress free before the temperature 
changes. Figures 5.50 and 5.51 show the results for the "original" specimen which does 
not contain TGO layer (i.e., h t = 0). The FGM top coat examined here is assumed to have 
linear material property variations, i.e., p x = pi = p$ = 1 in (5.6). It can be seen that 
FGM coating gives smaller crack opening (Figure 5.50) and strain energy release rate 
(Figure 5.51) compared with the homogeneous PSZ coating. The results which is 
consistent with the plane strain results given in Section 5.3.2.1 can, again, be explained by 
the lower thermal residual stresses and higher bending stiffness of the FGM coating. The 

* The material data in this table are the mean values computed from the temperature- 
dependent properties given in Table 1.1 while the data given in Table 5.4 are the same as 
the material properties at 22°C shown in Table 1.1. 
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results for the oxidized TBCs that consist of four layers are given in Figures 5.52-5.54. It 
is observed that the FGM coating with linearly varying thermomechanical properties gives 
smaller driving forces for interface crack propagation than homogeneous coating. The 
"rotated" stress intensity factors, as shown in Figure 5.54, indicate the stronger mode II 
contribution, which is also observed in the plane strain case. It is assumed in this study that 
the intrinsic residual stresses associated with TGO growth are negligible because of the 
expected stress relaxation due to the extended exposure to elevated temperatures. Note 
that because of the insertion of the TGO layer, the continuity of the material properties 
across the bond coat-top coat interface is lost and consequently the stress intensity factors 
are complex rather than real. Also note that as observed in experiments [27], the metallic 
particles in the graded coating would be oxidized after a sustained exposure to high 
temperature and, thus, a change of material property profile in the FGM is expected. 
Nevertheless, it was also been observed in the FGM coating system that the location of 
interface crack is right below the region which consists mainly oxidized metallic particles. 
Therefore a simplified layered system, i.e., FGM layer defined by (5.6) above the TGO 
layer is used in this study to approximate the more complicated real material system. 

Figures 5.55-5.57 show the effect of TGO layer on the crack opening and crack- 
extension driving forces for the TBC system with FGM (p 1 =p 2 =P 3 = lin (5.6)) top 
coat. It may be seen that the crack opens at a much smaller temperature drop when TGO 
layer exists. Also the crack driving forces for the case with TGO are much higher than that 
without TGO for the same temperature change. The huge differences in the results shown 
in Figures 5.55-5.57 induced by the thin layer of TGO may be attributed to the highest 
elastic modulus and lowest coefficient of thermal expansion among the materials in the 
TBC system (Table 5.5). This combination would increase the level of residual 
compressive stress significantly and thus gives higher strain energy release rate and larger 
crack opening. Therefore, it is very important to take the TGO layer into consideration in 
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evaluating the crack problems in TBCs when the surface temperature reaches beyond 
800°C. 

The strain energy release rates for the oxidized TBCs (Figure 5.48) with either 
homogeneous or FGM 0?i = P 2 = Pz — 1 in (5.6)) top coat are given as functions of 
normalized crack length in Figures 5.58 and 5.59. It is observed that under a uniform 


temperature change the strain energy release rate first increases rapidly but then gradually 
levels off as a/h c increases. An important application of the result of strain energy release 
rate as a function of crack length is to determine the service life of the coating under low 
cycle fatigue. It is well known that under cyclic loading such as the cyclic thermal stresses 
in the TBCs during service, subcritical crack propagation can be characterized as 


da 

dN 


= MAK), 


(5.7) 


where a and N denote the crack length and number of loading cycles, 
A K = iC max — AT,™ is the range of the stress intensity factor and the function / k is 
determined experimentally. For the mixed-moded interface crack problem discussed in this 
study, instead of (5.7), we may write 


da 

dN 


= fa(G), 


(5.8) 


where / G is determined experimentally and G may be expressed as G(a) (e.g., Figures 
5.58 and 5.59). By substituting the relation of G and a, exchanging dN and / G , and 
integrating both sides of (5.8), we may obtain crack length a as a function of number of 
loading cycles N. From this we may determine the number of fatigue cycles before a 
reaches the critical crack length at which brittle fracture occurs. 

Figures 5.60-5.65 show the normalized stresses for a homogeneous coating with TGO 
layer containing an interface crack of length a/h c = 30 (Figure 5.48) under 1000 °C 
uniform temperature drop. The normalized stresses for the specimen (Figure 5.48, 
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a/h c = 30) with FGM top coat (Pi = P 2 = P 3 = 1 in (5.6)) under 1000 °C uniform 
temperature drop are given in Figures 5.66-5.71. By examining the normalized stresses 
a rr , a zz , and a rz along the crack surface (0 < r < a) and the TGO-bond coat interface 
(a <r < R) for the cases of homogeneous top coat (Figures 5.60-5.62) and FGM top 
coat (Figures 5.66-5.68), it may be seen that, except for around the crack tip and the point 
of materials interface intersecting free surface, the residual bending stress a rT resulting 
from temperature drop is the dominant stress. It is also observed that the bottom surface 
of the debonded FGM coating is under higher compressive a rr than that for homogeneous 
top coat. The normalized stress a rT on the top surface of the homogeneous and FGM 
coatings are given in Figure 5.63 and 5.69, respectively. Compared to the stress a Tr on the 
top surface of the debonded FGM coating oy r for the homogeneous top coat appears to be 
higher and, thus, is more prone to the growth of surface crack. The normalized stress a rr 
along r = 0 and r = a for the TBC are given in Figure 5.64 and 5.65 for homogeneous 
coating and in Figure 5.70 and 5.71 for FGM coating. It may be seen that in both cases the 
TGO layer is under very high compression. 

The influence of material inhomogeneity on the crack opening at r = 0, strain energy 
release rate, and stress intensity factors are shown in Figures 5.72, 5.73, and 5.74, 
respectively, for the 4-layered TBC system described in Figure 5.48 (a/h c = 30). Here we 
assume that for the FGM top coat the inhomogeneity parameters for the elastic modulus, 
Poisson's ratio, and coefficient of thermal expansion have the same value, i.e.. 
Pi — P2 = Vz = V- Figures 5.75, 5.76, and 5.77 also plot the crack opening, strain energy 
release rate, and stress intensity factors as a function of inhomogeneity parameter p for 
several fixed values of temperature change. It can be seen from Figures 5.72-5.77 that the 
crack driving forces and opening are smaller for the metal-rich (p < 1) FGM than those 
for ceramic-rich (p > 1). Again, this is due to the better thermomechanical properties 
match, especially in the coefficient of thermal expansion, across the thickness direction in 
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the TBC system, and thus smaller residual stresses for the metal-rich case. Thus, by 
increasing the metal content in the FGM coating one may reduce the driving force for 
interface crack growth and consequently extend the service life of the TBC system. It is 
interesting to note that however, from the thermal insulation point of view, metal-rich 
FGM is less desirable because of its higher thermal conductivity. It is then important to 
take all relevant factors into consideration in order to obtain the optimal design of the 
FGM coating. 
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Figure 5.48 Specimen configuration for the second axisymmetric 
thermal buckling problem. 



Figure 5.49 The thickness variation of the Young's modulus in the 
FGM coating given in Figure 5.48, E\, and E c are the elastic 
modulus of bond coat and ceramics, respectively. 
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Figure 5.50 Crack opening displacements at the center of crack for homogeneous and 
FGM top coats without TGO under uniform temperature drop. To = 100 K. 



Figure 5.51 Strain energy release rate as a function of temperature drop for the 
homogeneous and FGM coatings without TGO, To = 100K. 
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Figure 5.52 Crack opening displacements at the center of crack for homogeneous and 
FGM top coats with TGO under uniform temperature drop. To = 100K. 



Figure 5.53 Strain energy release rate as a function of temperature drop 
for the homogeneous and FGM coatings with TGO, To = 100 K. 
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Figure 5.54 Normalized stress intensity factors for homogeneous and FGM coatings 
with TGO under uniform temperature drop, To = 100K, K 0 = E s a s TQ^Trh c , 

G 0 = (l-v s )KZ/E s ,K l + iKn = (K l + iK xl )K. 



Figure 5.55 Crack opening displacement at the center of crack for 
FGM top coat under uniform temperature drop. To = 100 K. 
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Figure 5.56 Strain energy release rate as a function of temperature 
drop for the FGM coating. To = 100 K. 



Figure 5.57 Normalized stress intensity factors for FGM coating under uniform 
temperature drop, To = 100 K, Kq — E s a s Toy/nh^, Go = (1 - u s )Kq/E s , 

Ki + iKn = (Ki + iKn)h«. 
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Figure 5.58 Strain energy release rate as a function of normalized crack 
length for the ho mog eneous coating with TGO layer, To = 100 K, 

K 0 = EsCtsToy/nh^, G 0 = (1 - v s )K$/E s . 



Figure 5.59 Strain energy release rate as a function of normalized 
crack length for the FGM coating with TGO layer, To = 100 K, 
Kq = E s a s Toy/nh~ c , G 0 = (1 - v s )K$/E s . 
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rla 

Figure 5.60 Variation of normalized 0Vr(r, zq ) along the crack surface (0 < r < a) 
and the BC-TGO interface (r > a) of the homogeneous TBC described in Figure 
5.48, z 0 /h c = 25.36, a/h c = 30, T 0 = 100 K. 



rla 

Figure 5.61 Variation of normalized o zz (r, zq ) along the BC-TGO interface 
of the homogeneous TBC described in Figure 5.48, zo/h c = 25.36, 
a/hc = 30, To = 100 K. 
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rla 

Figure 5.62 Variation of normalized a TZ (r, z 0 ) along the BC-TGO interface of the 
homogeneous TBC described in Figure 5.48, z 0 /h c - 25.36, a/h c = 30, T 0 = 100 K. 



Figure 5.63 Variation of normalized a TT (r, zq) along the coating surface of the 
homogeneous TBC described in Figure 5.48, z 0 /h c = 26.42, a/h c = 30, T 0 = 100K. 
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Figure 5.64 Variation of normalized a rr along r = 0 in the homogeneous 
TBC described in Figure 5.48, a/h c = 30, To = 100K. 



z/h c 

Figure 5.65 Variation of normalized o TT along r /a = 1 in the homogeneous 
TBC described in Figure 5.48, a/h c = 30, To = 100 K. 









rla 

Figure 5.66 Variation of normalized a TT (r, z 0 ) along the crack surface (0 < r < a) 
and the BC-TGO interface (r > a) of the FGM coating described in Figure 5.48, 
zo/h c = 25.36, a/h c = 30, To = 100 K. 



Figure 5.67 Variation of normalized a zz (r, z 0 ) along the BC-TGO interface 
of the FGM coating described in Figure 5.48, z 0 /h c = 25.36, a/h c = 30, 

To = 100 K. 
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Figure 5.68 Variation of normalized a TZ (r, zo) along the BC-TGO interface 
of the FGM coating described in Figure 5.48, zo/h c = 25.36, a/h c = 30, 

To = 100 K. 



Figure 5.69 Variation of normalized a TT (r, zq) along the surface of the FGM 
coating described in Figure 5.48, zo/h c = 26.42, a/h c = 30, To = 100 K. 





Figure 5.70 Variation of normalized cr rr along r = 0 in the FGM 
coating described in Figure 5.48, a/h c = 30, To = 100 K. 



z/h c 

Figure 5.71 Variation of normalized o TT along rfa = 1 in the FGM 
coating described in Figure 5.48, a/h c = 30, To = 100 K. 
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Figure 5.72 The influence of inhomogeneity parameter p on the crack opening 
displacement at the center of crack for FGM top coat with TGO under uniform 
temperature drop, T 0 = 100K, a/h c = 30. 



Figure 5.73 The influence of inhomogeneity parameter p on the strain energy release 
rate for FGM coating with TGO under uniform temperature drop. To = 100 K, 
a/h c — 30. 
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Figure 5.74 The influence of inhomogeneity parameter p on the stress intensity factors 
for the FGM coating with TGO under uniform temperature drop. To = 100 K, 
a /h c = 30, K 0 — E s a s Toy/ 7rh c , Gq = (1 — v s )Kq / E s , K\ + iTTn = (K\ + i-^n)^c € - 



ln(p) 

Figure 5.75 Crack opening displacement at the center of crack as a function of 
inhomogeneity parameter p for FGM top coat with TGO under uniform temperature 
drop. To = 100 K, a/h c = 30. 
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In (p) 

Figure 5.76 Strain energy release rate as a function of inhomogeneity parameter p 
for FGM top coat with TGO under uniform temperature drop, To — 100 K, 
a/h c = 30. 



In ip) 

Figure 5.77 Stress intensity factors as functions of inhomogeneity parameter p for 
FGM coating with TGO under uniform temperature drop, To = 100 K, a/h c = 30, 
Kq = E s a s Toy/irh c . Go = (1 — u & )Kq/E s , K\ + \Ko = {K\ + iKu)h 1 ^. 
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As discussed above and in Chapter 1, the oxidation process occurring at elevated 
temperature replaces the fabricated continuity with a distinct jump in material properties, 
introduces a weak interface, and therefore, accelerate the failure process of the TBC 
system. A more pragmatic approach to the design of TBC system for high temperature 
environment is to incorporate additional layers in the coating structure to satisfy the 
functionality required for more oxidation resistance. One example of such an approach 
would be to replace the original thermal insulation ceramics/metal grading by a multi¬ 
layered grading that consists of a thermal insulation ceramics/oxygen barrier ceramics 
FGM and an oxygen barrier ceramics/metal FGM. To study the blister and spallation 
problem for such TBC system, an axisymmetric specimen as shown in Figure 5.78 is 
considered. In this example, AI 2 O 3 is used as the oxygen barrier material and the same 
materials in previous example as described in Table 5.5 are used for the alternative TBC 
system. It is assumed that a thin layer of TGO would grow between the two FGM layers 
via the inward diffusion of oxygen through the PSZ/A1 2 0 3 FGM and the outward 
diffusion of cation toward the surface of the Al 2 0 3 /bond coat FGM. The dimensions of 
the specimen described in Figure 5.78 are given as: h c = 130 /tm, h g = 10/mi, h b = 50 
/an, h s = 3 mm, 2 R = 25.4 mm. A crack of radius a/h c = 30 is assumed to be along the 
TGO-Al 2 C> 3 /bond coat FGM interface. The temperature drop induced buckling problems 
are investigated for the "virgin" specimen (i.e., h t = 0 in Figure 5.78) and the specimen 
under extended exposure to high temperature environment ( h t = 5 /im). The 
thermomechanical properties of the PSZ/A1 2 0 3 and Al 2 0 3 /bond coat FGMs are defined as 
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respectively. In Equations (5.9a) and (5.9b) (E c , v c , a c ), (E t , v x , a t ) and (Eb, Vb, <*b) are 
the thermomechanical properties of PSZ, AI 2 O 3 , and the metallic bond coat, respectively, 
and are given in Table 5.5. Note that the power 4 in (5.9a) is used for a PSZ-rich top coat; 
and the power 1/4 in (5.9b) indicates a metal-rich FGM. The results are given in Figures 
5.79-5.81. For the purpose of comparison. Figures 5.79-5.81 also plot the results for a 4- 
layered material system with ceramics/metal FGM top coat described in Figure 5.48 with 
dimensions given as: h c = 130 /xm, h t = 5 fx m, hb = 50 /xm, h s = 3 mm, 2 R = 25.4 mm, 
a/h c — 30. The material inhomogeneity for the FGM top coat in the 4-layered system is 
defined in (5.6) with that pi = P 2 = P 3 = P = 4. It may be seen from these results that 
the crack opening at z = 0 and the strain energy release rate are much higher for the TBC 
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system with two graded layers than the original one with only FGM top coat. This is due 
to higher degree of mismatch in the coefficient of thermal expansion caused by higher 
content of aluminum oxide. Therefore, from the viewpoint of crack driving forces, one 
may expect that the multiple layered TBC system containing A1 2 0 3 would not provide 
better protection against coating blister and spallation. However, such coating system 
possesses other advantages against coating fracture. One such advantage is that the multi¬ 
layer arrangement eliminates the unfavorable material property discontinuity and 
consequently, removes the stress concentration and reduce the possibility of crack 
initiation at the locations where the interface intersects the free end [6], Also because of 
the higher density of the fabricated A1 2 0 3 oxygen barrier than that of thermally grown 
oxide, it provides better insulation against oxygen diffusion. This, in turn, reduces the 
growth speed of TGO and retards the weakening process of the "interface" (initiation and 
growth of microcracks) and, as a result, prolongs the service life of the coating [19]. It is 
then important to consider from both the viewpoints of increasing the material toughness 
and reducing the fracture driving forces when evaluating a coating system such as the 
TBC system with multiple graded layers described above. 
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Figure 5.78 Specimen configuration for the TBC with multiple 
graded layers under thermal buckling. 



A T/T 0 

Figure 5.79 Comparison of the crack opening displacement at the center of crack 
for the TBC systems described in Figures 5.48 and 5.78 under uniform temperature 
drop. To = 100 K. 
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A T/T 0 

Figure 5.80 Comparison of the strain energy release rates for the TBC systems 
described in Figures 5.48 and 5.78 under uniform temperature drop, To = 100 K. 
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Figure 5.81 Comparison of the stress intensity factors for the TBC systems 
described in Figures 5.48 and 5.78 under uniform temperature drop. To = 100 K, 


K 0 = E s a s T 0 ^h c , G 0 = (1 - v s )K$/E s , Ki + i K u = ( K , + iK u )h«. 
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5.4 The Influence of Initial Curvature 

Thermal barrier coatings are often used to protect structural components such as 
turbine blades in high performance engines and wing-tips for aerospace plane from high 
temperature environment. In these applications usually the coatings are applied to surfaces 
with curvature rather than flat. In this section we examine the effect of the initial curvature 
of the specimen on the blister of interface crack by considering a simplified model using 
the geometrically nonlinear finite element method. As an example, the plane strain problem 
of a ceramic-rich FGM coating bonded to a cylindrical metallic substrate containing an 
interface crack as shown in Figure 5.82 is considered. Thickness of the coating and 
substrate and the crack length are given as: h c = 130 [im, h s = 3.05 mm, and 2a = 5.2 
mm(a/h c = 20), respectively. The thermomechanical properties of the ceramic-rich 
FGM are defined as follows: 

E(z) = E C + (E m - < c +^ ) 4 , 


u(z) = v c + (y m - , R < r < h c + R, (5.10) 


a(z) = a c + (a m - a c ) (“ £ -^—-) > 


where the subscript c and m denote the ceramics and metal, respectively, and the material 
properties of the metal and ceramics are given in Table 5.6. The influence of curvature is 
studied by varying the radius of curvature of the cylindrical specimen. By taking advantage 
of the symmetry, only the right half of the cylinder is modeled. Since the purpose of the 
study is to examine the effect of the curvature, for simplicity, the bond coat and the oxide 
layers are not included in the analysis. 
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Table 5.6 Material properties of the cylindrical specimen described in Figure 5.82. 



Material 

E (GPa) 

V 

a (10' 6 °C' 1 ) 

Metal (Substrate) 

Ni-based Superalloy 

151.4 

0.26 

15.38 

Ceramics 

Zr0 2 -8wt%Y 2 0 3 (PSZ) 

11.54 

0.25 

10.91 


Figures 5.83-5.85 show the results for the specimen described in Figure 5.82 under 
uniform temperature drop. The constants used for normalization in these figures are 
To = 100K, Kq = E s a s Toyfnh c , Go = (1 - vl)Kl/E s . It is assumed that the cylinder 
is stress-free before the temperature change. The crack opening at x = 0 is shown in 
Figure 5.83 in which the case of R — infinity is approximated by a finite specimen as 
shown in Figure 5.21 with If a — 20. It may be seen from Figure 5.83 that, for the cases of 
finite radius of curvature R, instead of remaining closed until certain critical temperature 
change, the crack opens as soon as temperature starts to drop. Figure 5.83 also shows the 
results based on small deformation (linear) theory. As expected, the results for the linear 
model are simply straight lines in the plot of crack opening vs. temperature drop. By 
comparing the results based on linear theory with the one based on nonlinear theory, one 
may see that the results are similar for small radius of curvature; and that for large values 
of R the linear model predicts only the portion corresponding to initial temperature drop. 
The same observation can be made from Figure 5.84 in which the strain energy release 
rate is given as a function of temperature drop. It may also be seen that the results for the 
case of R = 2000 shown in Figures 5.83 and 5.84 behave basically the same way as the 
imperfection problems in buckling analysis. From this example we may conclude that for 
the problem of crack at a curved interface under in-plane compressive loading, the 
deformation is governed by nonlinear kinematic relationship when the curvature of the 
interface is very small; and that as the curvature increases (or R decreases), the 
deformation changes gradually from the one associated with buckling mechanism to the 

one governed by linear kinematic relationship. 
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Figure 5.85 shows the stress intensity factors as a function of temperature drop for 
several radii of curvature. It may be seen from this figure that the negative mode II stress 
intensity factor K\\ is the predominant stress intensity factor and increases (in magnitude) 
as the radius of curvature R decreases. Also observed in Figure 5.85 is that the mode I 
stress intensity factor decreases as R decreases. Furthermore, as the value of R drops 
below one hundred (results not shown here), K\ becomes negative and the crack surface 
contact near the crack tips is observed as soon as the temperature drops. To further 
examine the crack problems with such R values, a finite element procedure with contact 
algorithm is required. The crack problems with surface contact, however, will not be 
discussed in this study. 
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(a) The cross-section of the cylinder with an interface crack 


y 



(b) Enlarged crack region 

Figure 5.82 Configuration of the cylindrical specimen for the temperature loading 
problem, (a) the cross-section of the cylinder, (b) enlarged crack region, 2 a = R6. 
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Figure 5.83 Crack opening displacement at the center of crack for 
FGM coating under uniform temperature drop, To = 100K. 



0 5 10 

A T/T 0 

Figure 5.84 Strain energy release rate as a function of temperature 
drop for FGM coating, To = 100 K. 
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Chapter 6 

Conclusions and Future Work 


6.1 Conclusions 

In this study, the interface crack problem of TBC systems with functionally graded 
coatings under thermal or mechanical compressive stress is examined. The failure mode is 
characterized by the buckling instability induced coating blister which provides the driving 
force for the interface crack to extend and lead to the spallation of the coating. The 
instability problem for a FGM coating subjected to a far-field fixed-grip compression is 
investigated first by reducing the nonlinear governing equation to a linear eigenvalue 
problem using a perturbation method. The problem is solved under the mixed-boundary 
conditions by using the integral transformation technique. The analytical solution is then 
used to investigate the influence of coating inhomogeneity upon the instability load for the 
coating with either an interface crack or a highly weakened interface region and serves as 
a benchmark for the subsequent numerical study. 

The postbuckling problem is examined by using a geometrical nonlinear finite element 
method with enriched crack-tip elements for extracting the postbuckling fracture 
mechanics parameters such as strain energy release rate and stress intensity factors. The 
finite element procedure is used to investigate both plane strain and axisymmetric 
problems with various layered structure subjected to in-plane residual compression 
induced by uniform temperature drop. Solutions based on plate theory are also given to 
compare with the results based on the continuum theory. 

From the results obtained in this work, we may summarize the conclusions as follows: 
(1). For a FGM coating under in-plane fixed-grip compression as described in Figure 2.1, 
the instability strain is the highest when the coating has approximately the same modulus 
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as the substrate, and the value of the instability strain decreases for the case of stiffening 
( 7 h > 0) as well as softening ( 7 h < 0) FGM coatings. 

(2) . The instability strain for the graded coating with weakened interface zone is higher 
than that for the coating with a fully separated crack of the same size. The increase in 
instability strain is more significant for higher weakened interface length/coating thickness 
ratio. It is also observed that due to a relatively lower bending stiffness but the same 
bridging spring constant, the FGM coating with negative 7 h requires much higher 
compressive strain to induce instability compared with that for the coating with positive 
7 h. 

(3) . From the postbuckling analysis it is observed that the mode II stress intensity factor is 
negative and is the dominant stress intensity factor for the interface crack problem under 
in-plane thermal or mechanical compressive loads. From the viewpoint of crack driving 
forces, this means that the interface crack tends to kink into the coating rather than 
substrate and could result in the spallation of the coating. However, initiation of the crack 
growth requires that the crack driving force exceeds the mode mixity-dependent fracture 
toughness around the crack tip. Therefore, it is possible that the crack will extend along 
the rather weak (in toughness) interface until it reaches the free edge or until the mode II 
driving force is strong enough to cause the crack kinking. 

(4) . Compared with solutions based on continuum theory, the plate theory based analysis 
overestimates the instability strain and underestimates the strain energy release rate due to 
the boundary conditions with higher constraints. One may conclude that the plate 
approximation gives nonconservative results. The difference between results from the 
plate and continuum models lessens as the crack length/coating thickness ratio increases. 
Nevertheless, as observed in one of the examples, the postbuckling solution obtained by 
using plate model is limited to the first buckling mode, while under the continuum 
mechanics approach higher buckling modes are permitted. 
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(5) . Because of the better match in the thermomechanical properties in the FGM 
coating/metallic bond coat/substrate system compared with that for homogeneous ceramic 
coating/bond coat/substrate system, in the FGM coating the in-plane thermal residual 
compressive stresses is lower. Also the FGM coating possesses higher bending stiffness 
than pure ceramic coating. As a result the FGM coatings have certain advantages over 
homogeneous coatings from the viewpoint of crack driving force. 

(6) . Higher metal content in the FGM coating further improves the thermomechanical 
property match and increases the bending stiffness of the coating and consequently, further 
reduces the crack driving force. Nevertheless, higher metal content compromises the 
thermal insulation function of the coating. Therefore, the optimal design of the FGM 
coating would be to design the composition such that the coating provides the best 
fracture resistance while it satisfies the thermal insulation requirement. 

(7) . The TGO layer, with it's low coefficient of thermal expansion, is considered as the 
unfavorable composition in the TBC system because it increases the level of residual 
stresses and introduces a distinctive interface to the originally tailored FGM coating 
system. Besides, the oxide layer often served as the weak cleavage plane for the crack to 
initiate and grow. 

(8) . The problem of compression induced interface crack growth for FGM coatings on 
non-flat structural components is investigated. It is observed that for the case of small 
curvature the deformation of the coating is governed by the imperfection-buckling like 
mechanism and requires nonlinear analysis. On the other hand, for higher substrate 
curvature the problem may be characterized by linear small deformation theory. 

6.2 Remarks on Future Work 

The spallation of TBCs is basically a thermal-chemical-mechanical process. While the 
compression induced coating blister that leads to spallation has been investigated in this 
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work, there are still several ambiguous points in the whole picture of the spallation failure 
of the coating that need to be investigated. The first one would be, as briefly discussed in 
Chapter 1, the growth of TGO at the bond coat/top coat interface and its influence on the 
initiation of cracks along the interface and on the weakening of the fracture toughness of 
the interface. It has been observed that the spallation of the coating is strongly linked to 
the growth of oxide. The growth of TGO would result in the change of stress state near 
the oxide layer due to the volume change associated with the chemical reaction. Also the 
growth of TGO may introduce interface asperities at which cracks would initiate under 
local tensile stresses. Besides, the change in the constituents and microstructures of the 
TGO seems to have strong influence on the weakening of the fracture toughness of the 
bond coat/top coat interface. Therefore, the following problems are to be considered: (i) a 
model for the growth process of TGO and the associated residual stresses, (ii) a stress 
analysis and the associated crack problem for the interface asperities considering the effect 
of oxide growth, and (iii) a quantitative assessment of the interface fracture toughness as a 
function of oxide growth. 

One important problem to be investigated is to examine the influence of crack surface 
contact. It is observed that when the in-plane compressive load reaches a certain level in 
the postbuckling regime, crack surface contact occurs in front of the crack tip. It is also 
possible that the coating might deform into higher buckling mode which may cause crack 
surface contact. To correctly interpret the postbuckling results it is necessary to implement 
the contact algorithm and a special enriched crack-tip element which contains the correct 
asymptotic displacement fields considering crack surface contact. The crack surface 
contact also plays an important role in the more realistic three-dimensional coating blister 
problem, which can be an extension of this work. As observed in the delamination problem 
of composite laminates, a circular blister would gradually transform first to an elliptical 
shape and further into more complicated configurations [55,56,62] as the in-plane 
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compressive load increases. This feature, which cannot be detected by using a two- 
dimensional or axisymmetric analysis, is due to the crack surface contact near the crack 
front. Therefore, the contact algorithm and special enriched elements are required in the 
three-dimensional analysis to obtain the correct results. 

In this work we assume that the material behaves linearly elastically during the cool 
down period and that initially at high temperature the residual stresses are relaxed for the 
fracture characterization. In reality, because of the complexities in environmental and 
loading conditions during the processing and service of the TBC system, it is expected that 
the nonlinear material behavior would play an important role in the stress analysis and 
fracture characterization. In this work, as a first step, the fracture mode characterization 
provided an insight into the crack problem through the associated parametric study under 
the elastic assumption. However, to further examine the more realistic problems it is 
required to take into account nonlinear material properties. Those to be considered include 
the creep behavior of the materials during the period subjected to sustained exposure to 
high temperature, the sintering of ceramic top coat and the accompanying material 
property change, and the plastic deformation at the highly strained locations. 

Another problem to be considered as an extension of this work is to take into account 
the heat transfer in the thermal stress and crack problem. During the service stage, the 
temperature gradient and local hot spots in the material system are observed due to the 
heating and cooling boundary conditions and the local thermal insulation associated with 
cracks. The thermal stresses associated with the multidimensional temperature profile 
might be the driving force for the growth of interface crack. Also the non-uniform 
temperature change in the medium during the cool down process would result in different 
thermal residual stress profile than that for uniform temperature change. It is then 
important to consider the heat transfer in order to more accurately quantify the fracture 
driving forces. 
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Appendix A 

Finite Deformation Theory of Elasticity 


A.l Analysis of Strain 



Figure A. 1 Motion of a continuous medium. 


In characterizing the strains in the finite deformation theory, we employ the 
Lagrangian approach, in which the coordinates defining a point of the medium before 
deformation are used for locating the point during the subsequent deformation. Figure A.l 
shows at configuration t a deformed body t lR which occupies the region °R in the 
undeformed state. The coordinates of a typical material point p in the undeformed state 
are (°Xi, °x 2 , °x 3 ). The same material point occupies (*xi, *x 2 , f x 3 ) in the deformed 
state and we have 

Xj = Xj ~h zii, i — 1, 2, 3. (Al) 

Note that the left superscript is used to denote the configuration. To determine the strain 
in the Lagrangian configuration, first we define the deformation gradient [76] referred to 
the undeformed configuration as follows: 
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/ d*xi 

d*X! 

d*xi \ 

d°x\ 

d°x 2 

d°xs 

d*x 2 

d t x 2 

d t X 2 

d°xi 

d°x 2 

d°x 3 

d*X3 

d f xz 

d*X3 j 

\ d°x\ 

d°x 2 

d°X3 ) 


(A2) 


and define the spatial deformation gradient [76] which referred to the configuration t 
(Eulerian formulation) as 



of” 1 = {?*«} 


*0.13x3' 


(A3) 


Note that in the notation we used while a differentiation is involved, the left subscript 
indicates the configuration in which the coordinate is measured, for example. 


d ^ X' 

b x i,j = 0 n 1 . For an arbitrary infinitesimal material vector d°x at undeformed 


t nr . . — 


d°xj 


configuration and the corresponding vector d l x at deformed configuration we may write 

d l x = o F»d°x, d°x = o F~ l9 d l x. (A4) 


The strain tensor q E are defined such that it gives the change in the squared length of the 
material vector d°x as 


(d f s) 2 - (d°s) 2 = 2{d 0 xy 0 E»(d 0 x), 


where 


(Jsf = (d'xyid'x), (< d° s f = (d°xy(d°x). 

From (A4)-(A6) it may be seen that 


or 


t 1 f d^Xk dy_Xk 
° 6ij ~ 2 \d°Xi d°Xj 



i, j = 1, 2, 3, 


(A5) 


(A6) 


(A7) 


(AS) 


164 



where I is the unit tensor and <5^ is the Kronecker delta. By substituting (Al) into (A8), 
we may write the Lagrangian strain components as 


n (oUi,j ~ o Vj,i + o'U-k.i Q^k.j)- 


Note that (A9) is also called Green-Lagrange strains. Also note that if the displacement 
gradients are small, we may neglect the last term in (A9) and recover the components for 
small strains. 


A.2 Stress and Equations of Equilibrium 



Figure A.2 Force vectors for Piola-Kirchhoff stress definition. 


The stress measure to use with the Green-Lagrange strain tensor \E is the second 
Piola-Kirchhoff stress tensor q<t. The basic ideas of the second Piola-Kirchhoff stress 
tensor can be indicated in terms of the force vectors illustrated in Figure A.2. The second 
Piola-Kirchhoff stress tensor Icr is formulated such that it gives a force d°p related to the 
force dtp in the same way that a material vector d°x at °» is related by the deformation 
to the corresponding spatial vector d 1 x at f x. That is 


165 


(A10) 


( 0 n* t 0 <r)d°S = d°p = qjP 1 »d t p= t 0 F 1 • (n» t r)d t S, 

where V is the Cauchy stress tensor (defined in the Eulerian configuration). To determine 
the relation between Cauchy stress tensor and the second Piola-Kirchhoff stress tensor, 
first we have the relation between the deformed and original configurations given by [76] 

(All) 

d*V J ,, 
d°V ~ det{ ° n 


giving 

(A12) 

or 




t t t 

O x i,m O x j,n O^mn* 


(A13) 


Note that the second Piola-Kirchhoff stress tensor is work-conjugate with the Green- 
Lagrange strain tensor [76]. Also note that from the symmetry of V, it may be seen from 
(A12) that the second Piola-Kirchhoff stress tensor qCt is also symmetric. 

From (A10) we can also write the natural boundary condition in Lagrangian 
formulation as 

= (A14) 

or 

(o°i j + o a jko u i,k) ° n J = oTi> * = 1. 2,* 3, (A15) 

where qT is the applied traction vector. Note that the left subscript in qT, denotes the 
configuration with respect to which the quantity is measured. 
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From the static equilibrium of a deformed infinitesimal parallelepiped d*V, we may 
write the equilibrium equation for the Eulerian configuration as 

|r 2 + 'Pi = 0, i = 1,2,3, (A16) 

where t T ij is the component of Cauchy stress tensor and l P t is the Eulerian component of 
body force vector. To transform (A16) into Lagrangian formulation, first we consider the 
body force part in (A16). Assuming that originally the infinitesimal parallelepiped has 
volume d°V at undeformed configuration, we may write the relation between t P i and its 
equivalent body force component l Q Pi at undeformed configuration as follows: 

'Pid'V = iPid°V. {All) 

By substituting (A13) and (A17) into (A16) and using chain rule, we may write the 
equilibrium equation in terms of the second Piola-Kirchhoff stresses as 

<9 

q o ( 0%i,Tn 0%j,n Q^mn) t ^k,j T" oPi = 0- (A18) 

O Xh 

Providing that 

d ( d l Xj \ d°s fc = 

d°Xk \d°x n ) d t Xj 

3 (A19) 

O^j.n t't'kj = &kny 

Equation (A18) may be rewritten as 

oq ( Q&mn) Pi ~ 0. (A20) 

o x n 

The equilibrium equation (A20) for the Lagrangian formulation may further be simplified 
as 

-Q^r o'U'ijko&jkJ "f" -ft = ® = 1) 2, 3. (A21) 
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A.3 Stress-Strain Relations 

By assuming that the material is isotropic and generalizing the linear elastic relation for 
small deformation, we may obtain the elastic material description for finite deformation 
analysis as follows: 

0°ij = 0 Cijrs o ^rs dijCrer, i, j = 1,2, 3, (A22) 


where 


0 C,jrs — X6ij8 rs + (J,{8 iT 8j S "b 8{ S 8j T '), 

Ev E 

l “ (1 + i/)(1-2i/)’ M- 2(1+i/)’ 

C T = Y~2,v’ 6t = 


(A23) 


E, v, and a being the Young's modulus, Poisson's ratio, and thermal expansion 
coefficient, respectively, AT being the temperature change in the medium. 
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Appendix B 

Expressions of Various Functions in the Kernels 


B.l Expressions of A\, A 2 , C\, C 2 , 03 , and C 4 in Terms of F\ and F 2 

The unknown coefficients A\, A 2 , C 4 for the general solution in the stability 
analysis may be expressed in terms of F\ and F 2 , the Fourier transform of the density 
functions fi and f 2 , respectively, through a system of linear equations, (2.39)-(2.44), 
derived from the boundary and continuity conditions. These equations are solved in such a 
way that at first, two unknowns, C 3 and C 4 , are determined in terms of the other two 
unknowns, C\ and C 2 , from a pair of equations ((2.39), (2.40)) which contains these four 
unknowns only, then applying the relation obtained from the first step into another pair of 
homogeneous equations ((2.41), (2.42)) to determine C 4 in terms of the last two 

unknowns, A\ and A 2 , which appear in the boundary conditions ((2.45), (2.46)) that used 
as the integral equations. Compared with solving the six unknowns directly in terms of F\ 
and F 2 , it can be seen later that this procedure provides a better algebraic structure for 
studying the asymptotic behavior of the kernels. 

In order to simplify the asymptotic analysis which will be described in Appendix C, 
first the following quantities are defined 


n 


* 

1 


a± 

\ a 


n 2 — 

% — 


n 2 

a 


a 






8^0 / (8e 0 ) 2 

(k + 1) 2 +P Y (re + l) 4 


_ 8 reeo__o / ( 8 e 0 ) 2 

(re + 1) 2 P V( /£ + 1 ) 4 


8**o o I (8e 0 ) 2 

(re + l) 2+P V (« + 1 ) 4 



(Bl) 


(B2) 


(B3) 
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n 4 _ p£ 

" 4 -W“ 'T + \| 


£ 2 , 1 8 « e 0 a I ( 8 e o ) 2 {% — K\^ /T>/1N 

4 +1 -(^-^V(^TIF-(^TtM (B4) 


* 

"life 


TO* 


2ri fc + p(3 - «Q£ 


(k + l)n£ 2 + p(« + l)£n£ - (« - 1) + 8 (^j) € o 


•, fc = 1,..., 4, 


(B5) 


_ Ai _ 8 e 0 

1 Id V is + 1’ 


(B6) 


Ao 8 (k - l)ep 

A2 |a| V 1 (« + l) 2 ’ 


(B7) 


where rrij, rij and Xj are defined by (2.26), (2.27), and (2.29), respectively, and 

e-j2|. 

|a| 7 a 


(B8) 


Following the procedure described above, we may rewrite (2.39) and (2.40), 
respectively, as 




(B9) 


5>; - mty*M*c k = o. 


(BIO) 


k= 1 


Solving (B9) and (B10) we obtain 


C 3 = (2l|!—+ (g4- ~^ )e^-“=WSc 2 , (Bll) 

\fl3O4 — G463 / •— GLaO '* / 


0^4 — < 2463 . 


C, = ( ai ^~ 83 ? 1 )eW-”Bhl»/{ Cl + C ^3-a 3 »2 \ e (,;-„;)| # /; C2j ( g 12) 

\a3O4 — <2463 / ~ G4O3 / 


0364 — a 4 63 . 


where 
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3 — K 


■;CS —.i 

(B13) 

Vo. 

II 

1 

V 

(B14) 


Equations (2.41) and (2.42) may be rewritten, respectively, as 


, \ 4 r 

2Ai + ^1 + Aj 2 J A 2 + ^2 


3 — K 

K — 1 


I _ * * 

+ m k Ti k 




C k = 0, 


From (Bll), (B12), (B15), and (B16), C\ and C 2 may be expressed as follows: 

( 2^22+^21^12 


Ci=- 

c 2 = ( 


'&11&22 — ^12^21 


\ a , ( ^12^22 ~ 2 A 2 k\2 \ . 
) Al + \knk22 — k\ 2 k 2 \ / 2 ’ 


^21 fell +2fc2i \ ^ + / 2A2fen — ^12^21 \ ^ 
k\\k 22 — k\ 2 k 2 \ / \knk 22 — ^ 12^21 


where 


ii2 = -1 - a; 2 , 

h\ = Ai + —, 


(B15) 


/ 1 \ 4 

(AI + —JAi + 2 A 2 A 2 — 'y~' j ( n k ~ m l)Ck = 0 . (B16) 

1 k~ 1 


(B17) 

(B18) 


(B19) 

(B20) 


&11 = d \ + CL $ | 

( &461 ” (Xi &4 ' 
\a364 — 6^463 / 

J e K-"l)l 7 |feA +a4 

/ d\bs — <3361' 

Va 3 6 4 - a 4 fe 3 . 

^ e (nI-nJ)| 7 |fc/^ 

(B 21 ) 

k \2 = CL 2 + &3 1 

/ a ^2 — < 3-2 64 ' 
^ 0-3 — <^463 - 

+ 04 

/ a 2 b ^ — 0362' 

'<2364 — <2463 


(B 22 ) 

&21 = &1 + ^3 ( 

'a 4 6i — ai6 4 > 
>. 0364 — (Z463 ^ 

) e K-* 5 )l 7 |fc/£ + & 4 ( 

'^63 -a 3 6i^ 
' O3 64 — 0463 / 

| e (nI-nJ)| 7 |fe/^ 

(B 23 ) 

&22 = 62 + &3 ( 

'a 4 &2 - « 2 ^ 4 > 
^364 — (3463 > 

j e (" 2 -« 3 )l 7 |A/C + ft 4 | 

' a 2 bz — azb 2 > 
^364 — (34 63 J 

| e ( n 2 _ 7 l 4 )l 7 l*A < 

(B 24 ) 
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Now, substituting (Bll), (B12), (B17), and (B18) into (2.43) and (2.44), which may be 
rewritten as 

, 1 v 4 77» 




(B25) 


we obtain 


— (A\ + Ai — Ci — C 2 — Cz — C 4 ) = —, 

ia 


^22 \ F\ ( d\2 \ F 2 


dud,22 — d2\d\2 


)h-( 


d\\d22 ““ C?21 d\2 ' i OL 


)-• 

/ ia 


(B26) 


(B27) 


where 


Ao = 


( rf 2i > 


f rfl1 ) 

i Fs 

^11^22 — <^21 d\2 ' 

a 

~ 1 

\dud 22 ” d2\d\2 * 

1 . 
ia 


\dibA — dth'k J 


(B28) 


+ |m 2 + m 3 ^- 

+ mi 


0364 — 04 63 / 

+ m * (\ e (nI-»J)| 7 l*/«] ( ~ 2fc 22 -*2lfcl2 
^364 — dibs' _ ^& 11&22 “ & 12&21 

/04 62 - fl2&4 \-(n^-nDHfe/^ 


^ \ &3 64 — (X4 63 / 

/ Q 2&3 ~ g 3 & 2 \ c («;-wi)h|fc/g / jglfcll + 2 fe 2 l 
Va 3^4 — 0463 / _ \&n &22 — ^ 12^21 


^ _1_ 

+ AT’ 


(B29) 


\a 3 D 4 — < 2463 / 


ra* |" Q 1^3 - <*361 Vfnr-wilMfc/f ( j 12^22 - 2A 2 fci 2 

V (^3^4 d4 63 ' V Ati 1 /[loo — fci o fcn 


+ ra 2 


+ W»s( 


+ m4 


\a 364 — < 2463 / J V&n &22 ““ ^ 12^21 

o 4 fe 2 -a2&4 V^-n;)H/ t /£ 

fl3&4 — CI4&3 / 

,* 7 Q 2^3 - 03^2 \ ( n ;-n;)Mfc/f /2 A 2 fcn - ligfcgi 

— dth? / \ fci i &oo — Atio&oi 


^11^22 — &12^21 


+ a 2 , 


(B30) 
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d 2 1 = fl + 

L \ 0,304 —0,403/ 

| |" a l fc 3 - Q3&1 / - 2fc 22 - ^2lfel2 \ 

Va 3 6 4 — < 2463 / J '& 11&22 - ^ 12^21 ' 

+ L /04^0264 \ (n5-«J)Hfc/« 

\a 3 b 4 -a 4 b 3 ) 

+ f °2&3 ~ q 3fc2 \ c (n;-7i:)Hfc/£l / ^21^11 + 2fc 2 l \ 

\a 3 b 4 - a 4 b 3 ) J \kuk22 - ki 2 k 2 \ / 

- 1 , 


(B31) 


d 22 = [1 + 

\a3b4 — 0403/ 

+ / a l&3 ~ fl3&l \ c fa?-wnMft/f1 f ^12^22 ~ 2A|fci2 \ 

\a364 — 0,463/ J \knk22 — ^12^21 ' 

T / d 4 b 2 - (l2b 4 \ {n ;_nt)h\h/Z ( B32 ) 

^ L \a 3 b 4 -a 4 b 3 J 

_j_ / dobs ~ Q 3 ^ 2 \ e (n«-n^)| 7 [/t/gl f gap'll ~ ^12^21 N 

Va 3 6 4 - a 4 b 3 ) \ \kuk 22 - k 4 2k 2 \' 

- 1 . 


By substituting (B27) and (B28) into (B17) and (B18) and then into (Bll) and (B12), 
Ci,..., C 4 may be expressed explicitly in terms of F\ and F 2 . 


B.2 Expressions of Dijk in the Kernels (2.48)-(2.51) 

In terms of the quantities defined in Appendix B.l we may rewrite (2.45) and (2.46) as 
follows: 


lim ~— 

y—>~0 27T 


A- 


2 Ait* 


iA2& >i2yS jicxe iax da = 0, — a < x < a, 


(B33) 


ton ^ f (hiA^v + 2\* 2 A 2 e X2y ) \ot\e iax da = 0, 


a < x < a. 


(B34) 


By substituting (B27) and (B28) into (B33) and (B34) the functions Dijk, i, j, k = 1, 2, in 
the kernels (2.48)-(2.51) are given by 
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(B35) 


An (a) 

Ai2(a) 

Dm(oc) 
A 22(a) 
Aii (a) = 

Ai2(a) 
A21 (a) 
A22(a) = 


2(^22 


d\ 1 

^22 — d2idn ’ 



dn 

d22 ~ d2idi2 ’ 


2di2 

dn 

d>22 ~ $21 d\2 ? 


li2dn 

dn 

&22 ~ <&2\d\2 5 


h\d 2 2 

~ d 

11^22 — d 2 idn 


2\* 2 d2i 

dn 

0?22 — d2ldi2 ’ 


hidn 

dn 

d 22 — d2idi2 ’ 


2A ld n 


dnd 2 2 — d2\d\2 


(B36) 

(B37) 

(B38) 

(B39) 

(B40) 

(B41) 

(B42) 


Note that from (B1)-(B5), (B13), (B14), (B21)-(B24), and (B29)-(B32) it may be 
seen that 


nj(a) = n^(-a), 

n* 3 (a) = n* 4 (- a), 

(B43) 

(a) =m* 2 (-a), 

m* 3 (a) = ml(-a), 

(B44) 

ai(a) = a 2 (-a), 

03(a) = a 4 ( - a), 

(B45) 

's' 

1 

<N 

II 

rO 

63(a) = 6 4 ( - a), 

(B46) 

^11 (a) = ki 2 ( — a), 

&2i (a) = fc 22 (-a), 

(B47) 

dij{oi) = dij( — 

a), i, j = 1, 2. 

(B48) 


Thus, from (B35)-(B42) we may see that 
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Dijk (ci) — Dijk{_ qj)> ol < oo, i , jy k — 1, 2. (B49) 


Furthermore, it may be seen that when |a/ 7 | < [(3 - /c)(l + /t) 3 ] 1 / 2 /8e 0 and eo < 0.5, 
n* (j = 1, ...,4) are complex. From (B1)-(B5), (B13), (B14), (B21)-(B24), and (B29)- 
(B32) we have 

n*j(a) = nf(-a), m*(a) = mj(-a), 

aj(a) = aj(~ a), bj(a ) = - a), i, j = 1, 2, (B50) 

M a ) = hj{ ~ <*)» dij(a) = dij( - a), 


when [or/'y| < [(3 — /c)(l + tf) 3 ] 1 / 2 /8eo- By applying (B50) to (B35)-(B42), we obtain 

[(3- K )(i + K y] 1/2 


Dijk (pc') Dijk( Qi), 


a 


7 


< 


8e 0 


i, j, k = 1, 2. (B51) 


It may be seen from (B51) that the real part of Ajjt's are even in a, which agrees with 
(B49). However, according to (B51) the imaginary part of D^k s are odd functions of a 
when |a/7| < [(3 — k)(1 + /c) 3 ] 1 / 2 /8eo, but (B49) gives that lm(Dijk) is even in a. 
Therefore, we may conclude that D^'s are real functions and are even in a. 
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Appendix C 

Asymptotic Expressions for the Integrands in the Kernels 


To asymptotically expand the integrands Bij, i, j = 1, 2, (Equation (3.10)) in the 
kernels for large 77 , we follow the approach for determine the closed form ((B35)-(B42)) 
of these integrands described in Appendix B. Beginning with the characteristic roots (Bl)- 
(B4), we determine the asymptotic expansions, substitute them into the expressions for the 
successive algebraic functions, and repeat the asymptotic expansion-substitution process 
until the asymptotic series for (B35)-(B42) are obtained. 

In order to simplify the asymptotic analysis, we first carry out the first-term asymptotic 
approximation for some functions. As a approaches infinity, from (3.5) and (B 8 ) we have 


V ->oo, £-0, (3 = 1. (Cl) 

Thus, from (B1)-(B5), (B13), and (B14), we have, as a or 77 goes to infinity, 

n;=-A$ + O(0, (C 2 ) 

n$=-AI+O(0, (C3) 

713 = A; + 0(£), (C4) 

»4=A; + O(0, (C5) 

m\ = \* 2 + 0(0, (C6) 

™2 = 4 + 0(0, (07) 

A i 

m* 3 = - A 2 + 0(0, (C8) 

«•; = - i+o«). (C9> 

A i 
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(CIO) 


a, = -1 - a; 2 + 0(0, 


a 2 — -2 + 0(0, 

(C11) 

a 3 = - 1 - Aj 2 + 0(0, 

(Cl 2) 

a 4 = ~ 2 + 0(0, 

(Cl 3) 

bi= — 2X1 + 0(0, 

(Cl 4) 

^2 = ~ ^1 _ ^7 + 0(0, 

(Cl 5) 

63 = 2 A 2 + 0 ( 0 , 

(Cl 6) 

6, = a; + i+ 0 ( 0 , 

(Cl 7) 

where A* and A£ are given in (B 6 ) and (B7), respectively. It may be seen that for 

1 < k < 3 and eo < 0.228, we have a+bj — afii = 0(1) (i, j = 1 , 2, 3, 4). Also it can be 

seen that the leading term for n\ 2 — 713 4 is negative and of order one. Therefore, from 

(B21)-(B24), we may write hj as 

hi = a x + 0(0), (Cl 8 ) 

&12 = 0-2 + 0(0), 

(Cl 9) 

hi = h + 0(0), 

(C20) 

&22 = 62 + 0 ( 0 ). 

(C21) 


Furthermore, from (B29)-(B32), the asymptotic expression for d tJ ( i , j = 1, 2) may be 
expressed as follows: 
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It is clear to see that instead of performing the asymptotic expansion for the extremely 
complicated (B29)-(B32), we may simply expand (C22)-(C25) to obtain the asymptotic 
expressions for dij. 

In this analysis, the 11-term asymptotic expression for A# and Bij, defined by (B35)- 
(B42) and (3.10), respectively, for 77 —> 00 or £ —*• 0 are desired. The procedure used for 
the analysis is, first, obtaining the 11 -term asymptotic expansions for n\ and n^, 
substituting them into (B5) to compute the 11 -term asymptotic expansions for m\ and 
m* 2 , obtaining the 11 -term asymptotic expressions for 01 , 02 , 61 , and b 2 by substituting the 
expansions of n\, n £, m\, and m* 2 into (B13) and (B14), then substituting the results into 
(C22)-(C25) to determine the 11-term asymptotic expansions for d^, and finally 
substituting these series into (B35)-(B42) and then (3.10) to obtain the 10-term 
asymptotic expressions for A# and Bij. It may be seen that the asymptotic analysis 
described above involves only simple algebraic manipulation of polynomials once the 
asymptotic expansions for n\ and n\ are obtained. 

The asymptotic analysis is performed by using symbolic manipulating software 
MAPLE. Because of the lengthy expression for the coefficients for higher order terms, 
only the expression for first 2 terms are given for A#' s and Bij s, which are given by 


A 11 


2A; 3A^ 2 + A^(l + AA* 2 -At 2 )-(l + At 2 ) 2 /p\ 

1-Ai 2 2Ag(Ai + A 2 XI — AJ 2 ) 2 VT/i 



(C26) 
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(1 + Af ) 2 
2 A^( 1 -Af) 

(A ^ 5 + Aj 4 A£ + 4Af + 3Aj - 8 AjAS 2 - A^)(l + Aj 2 ) /p 
8 Af(At + A^)(l-At 2 ) 2 \ V 

+ 0(~2 )> 

VTT / 


(C27) 


77i?i = — 


_ i + a ; 2 , 3A; - a; - a; 2 a; - a ; 3 /P \ , /i 


1-At 4A;(AJ + A5)(1 - At 2 ) V 77 . 


(C28) 


-D 122 = 


1 + At 2 (1 + At 2 )(l - AtA 2 ) /£ 

1 - At 2 4AtA£(At+A^)(l-At 2 )V77 


- +0 


(C29) 


Dm = — 


1 +Af 
1 - At 2 


(At 4 + Aj 3 A$ + 2At 2 - Aj 2 A $ 2 - At A?; - 3A $ 2 + 1)(1 + At 2 ) (p 


'(?)• 


4A;a;(a; + A5)(i-a; 2 )2 


0 (C30) 


1 + At 2 
1 - At 2 


_ (^i 5 + At 4 A^ + 4At 3 + 3At - 8AtA2 2 - A|) (P\ , n /J_\ 

4A^(AI+A^)(l-At 2 ) 2 V V'’ 


(C31) 


Dri 1 = — 


(1 + At 2 ) 2 
2At(l - At 2 ) 

(3A% — At — At 2 Ag — At 3 )(l + AJ 2 ) (p 


8 At 2 (At -h A^)(l - At 2 ) 


0 +O 0)' 


D,„ = __ (t 

i-a ; 2 2A;(a; + a;)(i-a; 2 )vii 


(C32) 




Bij - B ifl f ) + B ij2 ( 2 ) + 0 ( 0 ), j - 1, 2, 


(C34) 


where 
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(C35) 


Bm = p(Xf + XfX* 2 + 5Af - 3A? 4 A£ + 7Af - 12AJ 3 X? + 4Af A^ 3 


- 9Af A^ + ZXl - 

+ a^)(i - a; 2 ) 2 


5AS + 12A 


5 s )/[i 


8 A^(AJ 


[ 4 A 1 A 2 - (1 + At 2 ) 2 ] (A? 2 - 2 A 1 A 2 + l)(Af - 2Ag 2 + 1) 

112 8AjA5 3 (A; + x* 2 )(i - Aj 2 ) 3 

= p(At 2 + AfA^ 2 - 3A$ 2 + 1) 

121 4A^(AI+A^)(l-At 2 ) ’ 

[4A^-(l + At 2 ) 2 ](Af-2A^ + l) 

122 16XfXf(l-Xf) 2 


= p(Aj 2 + A; 2 A$ 2 -3A$ 2 + i) 
l_ 4A;A5(A; + ^)(1-AI 2 ) ’ 

[4A;A£ - (1 + Xf ) 2 ] (Aj 2 - 2A*As + 1) 

i6a; 2 ^ 2 (i-a; 2 ) 2 ’ 


_ p(A; 8 + A|^A| + 2Af - 6A? 2 A% + 5A? - 3 A 2 ) 
8Af(Aj+A" 2 )(l-Al 2 ) 

_ [4a;^-(i + a; 2 ) 2 ](a; 2 -2a;^ + i) 

222 8Aj s A5(A; + a^)(i-a^ 2 ) 2 


(C36) 

(C37) 

(C38) 

(C39) 

(C40) 

(C41) 

(C42) 


Note that the relation B\ 2 k = ~ B 2 \k, k = 1, 2,..., is observed in the asymptotic analysis. 
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Appendix D 

Formula for Evaluating the Asymptotic Kernels 


It may be seen that the integrals involving the asymptotic series in the Fredholm 
kernels (3.16), (3.17), (3.22), and (3.23) have the following two general forms: 


C k (u) = f 

JA 

S k {u) = 


COSLJTj 

r)& 


a n 


sin Ljq 

rnk 


dv, 

k = 1,2,..., 

(Dl) 

dih 

k = 1,2,.... 

(D2) 


Ja T 

To determine Ck and S*, first we make a change of variable as follows: 

t = M 7], 

and Ck and Sk can be rewritten as 

CiM = r -0-dt, 

J\w\A t 


_ . . |u )\ k f 00 sinf 

= / .jj. dt. 

u JUA t k 


Now consider the integrals in (D4) and (D5). For the case of k = 1 we have 

r cos t 


-dt = — Ci(|u;|A), 


r 

L * JO 


-dt — Si(|a;| A) 


= 5 - Si(MA), 


where 


. . , \ f cost — 1 , 

Ci(x) = 70 + ln(x) + / --- dt, x > 0, 

Jo t 
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Integrating by parts, we may express the integrals in (D4) and (D5) for k > 1 as follows: 


/ 

Jw 


\A 


COS t cosf 

= “ (k - 


1 f°° sinf 

\w\A k ~ 1 J\w\At k ~ l 

A) f°° sin t 

* -1 L(M^)*- 1 L\At k - 1 ' 


dt 


t—\uj\A 

cos(|o;|A) r °° 


<\uj\A 1 

-df] 


(DIO) 


sin t , 

/ -jrdt = - 
J\u>\A t 


sinf 


(k — l)t fc ' 


-l 


+ 




k 


l i* 00 

1 <-/ |£ jj | j 4 


cosf 

f fc_1 


dt 


sin(|o;|A) 

fe-iL(MA )*- 1 


r < ^L 

J\u>\A t k ~ l 


dt 


(Dll) 


By substituting (D6)-(D11) into (D4) and (D5), C*' s and S k ' s may be expressed as 
follows: 


C\((jS) = -Ci(|w|A), 


Si{u) - 


M 

u 


7T 


- Si (M^) 


Ci(w ) = r l_f£^Md) +wSl . lH 


, k > 1, 


S k {u) = 


As — 1 


u\ ^ sin(|o;|yl) 
A* -1 


a 


+ uCk-\ (w) 


k> 1, 


(D12) 


(D13) 


(D14) 


(D15) 



Appendix E 

Matrices and Vectors Used in 
the Two-dimensional Element Formulation 



Figure E.l Two-dimensional 12-noded isoparametric element. 

For the derivation of the matrices and vectors, we consider a typical two-dimensional 
element at time t as shown in Figure E.l. For isoparametric element formulation, we write 


the element coordinates and displacements as 


M 

M 

a x 2 = Y i N i °xl 

k =1 

k =1 

M 

M 

k =1 

k =1 


(El) 

(E2) 


where M is the number of nodes in the element, Nk is the interpolation function, °Xi and 
°#2 are the global coordinates of the nodal points at time 0, and t u\ and are the 
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displacements of the nodal points at time t. Similarly, the incremental displacements u\ 
and U 2 in terms of incremental nodal displacements u\ and u 2 are given by 

M M 

«i = J2 Nk M i 5 u * = ^N k u k 2 . (E3) 

*=1 k=l 

Note that the superscript ( n ) which refers to the values of displacements, strains, and 
stresses computed in the n-th iteration for certain increment is omitted in this Appendix. 
However, it is important to note that the values computed from (n - l)-th iteration are 
used for evaluating the matrices and vectors with a superscript (n). 

First we consider the axisymmetric problem in which the coordinate (°xi, °x 2 ) refers 
to (' r,z ). To determine the linear strain-displacement transformation matrix 0 the 
following incremental linear strain and incremental displacement vectors are defined 

oe (n) = |oeii, 0^22) 2oei2, 0633 }, (E4) 




where 
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(E7) 


t 

0 



/ 


V 


(1 + ln)oNi t i 

I 12 0^1,2 

(1 + Jll)oNl,2 + ^12 0-^1,1 

JVi 

(l + * 33 ) 5 ^ 


hi oN 1A 
(1 + ^ 22 ) 0 ^ 1,2 
(1 + ^22)0^1,1 + ^21 oiV 1(2 

0 


(1 + /n)o-Wi,i 
J 12 0^2 

(1 + hi)oNi^ +112 0-^,1 

(1+w % 


^21 oiV M 

(1 + 

(1 + h^oNi,! + hi oNi } 2 
0 


(1 + Iu)oNm,i 
ll 2 oNm ,2 

(1 + Iu)oNm,2 + Ii2oNm,1 

(1+w£" 


hi oNm,i 
(1 + 122 )oNm,2 
(1 + h2)^M,l + hi oNm,2 

0 


J Xi 


\ 


J 


1 < i < M, 


in which °xj denotes the °xi value of the Gaussian point where qB^ is evaluated, 0 Nkj 
is given by 


0 Nkj 


dN k 


J d°Xj' 


(E8) 


and kj are defined as 


M 


M 


E 


hj ^ ] o Nkj rtj, i, j 1)2, I 33 


fc =1 


fc=i 


°a;i 


(E9) 


Similarly, the nonlinear strain-displacement transformation matrix Bnl related to oWij in 
(4.36) is given by 
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(E10) 


Bnl = 


/oM,, 

0 

0 ^ 2,1 

0 

oN m ,i 

0 \ 

0 -Wl,2 

0 

0 ^2,2 

0 

()Nm ,2 

0 

0 

o N lfl 

0 

0 -^2,1 

0 

o Nm,i 

0 

oN\ ,2 

0 

0 -^2,2 ... 

0 

oNm ,2 


0 

N 2 _ 

0 

N m 

° J 

V °^i 

°Xi 

°Xi 


The stress matrix corresponding to B^l in (4.43) and the stress vector corresponding to 
in (4.45) may then be expressed, respectively, as 


ls in) = 


( 0^11 0^12 
0 a 12 0 a 22 
0 
0 

V o 


0 

0 


0 

0 


0 

0 

0 


0^11 0^12 
0^12 0 a 22 


0 \ 
0 
0 
0 


o o ()a 3 3} 


(Ell) 


i o( n )^_ Jt _ t t t \ 

0* — | 0^11) 0^22) 0^125 0^33 J- 


(E12) 


For the linear elastic, isotropic, and inhomogeneous medium investigated in this study, 
we may write the constitutive equation for the asymmetric problem as 

o'S'= C'o'e + q'S't, (E13) 


where OS' is given in (El2) and 


C = 


E 


(\ — v v 0 v \ 

v 1 — v 0 v 


(1 + j / )(1 — 2 v) 


1 — 2 v 

0 0 —— 0 


V 


V V 


0 l-v 


) 


l? 1 — |o e ll) 0 e 22) 2oCi2, o € 33 


t c»T EqT 

° St - ~ 0 T 2 u) 


{ 1, 1,0, l}, 


(E14) 


(E15) 
(El 6) 


T being the temperature change between "time t" and "time 0". Similarly, we may write 
the incremental thermal stress 0 St as 
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0 



Ea(AT) 

(1 - 2u) 


{l, l,0,l}, 


(El 7) 


AT being the temperature change between "time t + At" and "time t". 

For plane strain problem in which (°xi, °X 2 ) refers to (x,y), it may be seen that the 
matrices and vectors are very similar to the ones used for axisymmetric problems and only 
minor modifications on (E4)-(E17) are needed: for 0 e^ n \ qS^\ q€, qSt, and 0 St, we 
remove the last entry in the vectors in (E4), (E12), (E15), (E16), and (E17), respectively; 
for IB^ and Bjyjji the last rows in (E7) and (E10), respectively, are removed, for q^ 
and C, the last columns and last rows in (El 1) and (E14) are removed. 

The matrices and vectors used for the plane stress problem are the same as the ones 
used for the plane strain problem except the following constitutive law related matrix and 
vectors: 


C = 


1 — v 2 


(l v 

V 1 
0 0 


0 \ 

0 

l-v 

2 ) 


(plane stress), 


(El 8) 




3 t _ 
o^r — 


Hums 1 1 0 ) 

(1 - v) \ ’ ’ J 


(plane stress), 


(plane stress). 


(E19) 

(E20) 
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Appendix F 

Functions and Matrices Used in the Enriched Finite Elements 


The functions used in the asymptotic displacement expressions (4.51) for an interface 
crack between two distinct isotropic materials (Figure 4.2) are given by 

fij = F\j cos <f> - Gij sin (f>, g X j = F^j cos (j> - G^j sin (f>, 


foj — Fij sin (f> + G\j cos </>, g^j — Foj sin 0 + G 2 j cos <j>, 


(FI) 


where 


Fij = 


4 [ij cosh 


»sh ( 7 re) C 


+ 26 sin 6 sin 0), 


F2j= ~ 


Afij cosh ( 7 re) y 27 t 


— (.Ej - 26 sin 6 cos 0), 


Gl ’ = - j cosh(Ji )£ (E > + 2{sin#cos0 )- 

h^T)]fk {Di - 2Ssinesine) ’ 


G<ij = — 


4 jij cosh 


27T h 2 // \/i 2 


5 i = 6 2 = e (7r+(?)£ , 


(F2) 


(F3) 

(F4) 


© = elnr + 


6 d 6 6 

Dj = /?7jcos - + (3’^’ sin Ej = P'^cos - - sin 


(F5) 


(F6) 


a 0.5cos(elnr) + esin(elnr) a , 0.5sin(elnr) — ecos(elnr) /Ty7N 

^ ~ 025+^ ’ m 
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Bi(M+1) Bi(M+2) 

B 2 (M+1) B 2 (M+2) 

Bz{M+\) B 3 (M+2) 
B4(M+ 1) -S 4 (M+2) 
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B : 


3(JW+2) 




d°x 2 d 0 ^ 


_ y- <9iV fc fc &92j_ _ &N_k k 

dO Xn yij ao„ 02j 

k=l 


d°xi 


d°xi A] 
k =1 a Xl > 

M 


(FI 6) 




(FI 7) 



(FI 8) 


(F19) 


Note that the matrices (F9) and (F10) are given for axisymmetric problem. The linear and 
nonlinear transformation matrices for plane problem may be obtained by simply removing 
the last row in (F9) and (F10). The zeroing function Z has a value of unity for the 
enriched crack-tip elements; and for transition elements Z varies from one at the edges in 
contact with crack-tip elements to zero at the boundaries with regular nonlinear elements. 
A linear zeroing function Z in terms of local coordinate (r, s) is given by 


Z(r, s ) 


-(1 ± r)(l ± s), for zeroing from comer node, 

< -(1 ± r), for zeroing from element edge, 

- (1 ± s), for zeroing from element edge. 

2 


(F20) 
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Appendix G 

Test Examples of the Finite Element Implementation 


G.l The Elastica Problem 



Figure G.l (a) The elastica problem, (b) Finite element model 
used for simulating the elastica problem. 


The problem of large deflection of a buckled column (the elastica) [87] described in 
Figure G.l (a) is used to test the finite deformation implementation in the finite element 
code. In the elastica [87] it is assumed that the cantilever column under consideration is 
inextensible and the load P remains in the direction parallel to y-axis. The displacements 
at point A is given by [87] 
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(Gl) 


2 2 /? 
u„ = -£(/})-2(, v. = f, for P> Per, 

where 

k = yfp/EI., P = sin (a/2), kl = K(J3), P cr = 

(G2) 

/ >7r /2 f */2 _ 

K( 0 )= , ^ . , E(t3)= VTTpiftdt. 

Jo i/l-0 2 sin 2 4> Jo 

Per being the Euler buckling load, K{{3) and E({3) being known as the complete elliptic 
integral of the first and second kind, respectively. Note that the solution (Gl) is given in 
such a way that the unknown displacements {u a , v a ) and the input load P are expressed 
in terms of the deflection angle a at point A. 

Using finite element program, the model in Figure G.l(b) is used to study the elastica 
problem. To simulate the slender column the relative dimensions are used: l/h = 50, 
b/h = 1, where l, b, and h are the length, thickness, and width of the bar. By gradually 
increasing the uniform compressive pressure p ( = P/bh ), the displacements at point A 
are obtained as a function of P from the incremental-iterative finite element procedure 
described in Chapter 4. Intuitively it may be seen that under axial compressive loading 
condition as shown in Figure G.l(b) (Q = 0), the column would simply "shrink" rather 
than buckle no matter how large the compressive pressure is. Therefore, initially a 
transverse force Q acting as a disturbance to introduce transverse displacement of point A 
is applied to point A; a reverse loading - Q is applied at a later stage after the main 
loading P exceeded the critical load P a . By doing this we can obtain an equilibrium 
configuration for the column at certain value of P > P a on the preferred branch 
(nontrivial solution) in the displacement vs. loading bifurcation diagram. We may then 
start from this equilibrium configuration and gradually either increase P or decrease P to 
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recover the displacement/loading relation for the elastica problem. This solution procedure 
can be seen from Figure G.2 which gives the x- displacement of point A vs. loading P. In 
this particular case we choose the loading parameters P* = 1.063P cr and Q* = P* j 20. 
The axial and transverse loadings P, Q are incrementally applied in such a way that at first 
(i) increase P and Q from zero to P*/2 and Q*, respectively, then (ii) P — P*/2 —* P*, 
Q = Q* —>► 0, after that (iii) increase or decrease P from P* and Q = 0. The stages (i) 
and (ii) are shown in Figure G.2 as dotted and dashed line, respectively. 



Figure G.2 Transverse displacement of point A for the problem described in Figure 
G.l(b), the applied loadings: ®P = 0-> P* / 2, Q = 0 Q\ (ii) P = P*/ 2 — P*, 
Q = Q* ^ 0, (iii) P = P* -+ 3P CT or P = P* -> 0, Q* = 0 (P = bhp, 

P* = 1.063P ct , Q* = P*/20). 
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Figure G.3 The deformed shape of slender bar at P = bhp = 9 P c 



Figure G.4 Displacement components of point A in the elastica problem 
described in Figure G. 1. 





The deflected shape of the column at P/P a = 9 is shown in Figure G.3. Figure G.4 
shows the displacements of point A as a function of compressive load P from both 
analytical and finite element methods. As can be seen in Figure G.4, the results match 
quite well for 1 < P/P a < 1.5 and the difference between analytical and finite element 
results is less than 7% difference for 1.5 < P/P cv < 9.5. This difference may be attributed 
to the assumption used in the analytical formulation that the column is inextensible, and as 
expected, the difference becomes more significant at higher load. 

G.2 Interface Crack in FGM coating/Homogeneous Substrate 
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Figure G.5 Functionally graded coating bonded to homogeneous substrate with 
an interface crack subjected to uniform tensile stress normal to the surfaces, shear 
modulus of the coating: /z 2 (y) = //iexp (7 y). 


To test the inhomogeneous and enriched element implementation, we consider the 
problem of a FGM coating bonded to a homogeneous substrate containing an interface 


196 



crack subjected to uniform normal surface loading as described in Figure G.5 [72] where 
the Poisson's ratio v = 0.3 and the relative dimensions are given as: h\ = 4a, h 2 = a. It is 
assumed that the strain measures follow the small deformation theory and that the material 
is linear elastic. The strain energy release rates and stress intensity factors for various 
inhomogeneous parameters obtained analytically from singular integral equations [72] and 
computed by using inhomogeneous enriched finite elements are shown in Table D.l. Note 
that in the finite element analysis the specimen is finite in x direction and that to eliminate 
the end effect the length / = 40a is used. Again, excellent agreement between the 
analytical and finite element results are observed. 

From the results of the test examples in G.l and G.2, it is safe to assume that the 
inhomogeneous enriched finite element code with large displacement formulation is 
working correctly and that the finite element solutions for problems with unknown 
analytical solutions can be considered to be accurate results. 
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Table G.l Normalized stress intensity factors and strain energy release rate 
for the problem described in Figure G.5. 


7 

G/Gq 

h (a)/, 

M 1 '* 

k 2 (a)/ J 

P0<2 1/2 


Analytical 

FEA 

Analytical 

FEA 



3.0 

1.132 

1.136 

1.062 

1.064 

0.068 

0.067 

2.5 

1.244 

HBEB 

1.115 

1.117 

0.039 

0.038 

2.0 

1.382 

wmm 

1.175 

1.178 

0.006 

0.005 

1.5 

HESsH 

1.557 

1.245 

1.247 

KZS29I 



1 

1.766 

1.342 

1.327 

-0.076 

-0.076 

0.75 

1.880 

1.889 

1.367 

1.371 

-0.100 

-0.101 

0.5 

2.014 

2.024 

1.414 

1.417 



0.25 

2.164 

2.175 

1.463 

1.467 

-0.153 

-0.154 

0.1 

2.262 

2.273 

1.494 

1.498 

muiH 

-0.171 

-0.01 

2.337 

2.349 

1.518 

1.521 


-0.185 

-0.1 

2.401 

2.414 



-0.195 

-0.196 

-0.25 

2.514 

2.527 

1.571 


-0.214 

-0.215 

-0.5 

2.717 

2.732 

1.630 


-0.247 

-0.248 

-0.75 

2.942 

2.958 

1.692 

1.696 

-0.282 

-0.283 

-1.0 

3.191 

3.208 

1.758 

1.762 

-0.319 

-0.320 



3.483 

1.826 

1.831 

-0.358 

-0.360 

-1.5 

3.765 

3.786 

1.899 

mm 

-0.400 

-0.402 

-2.0 

4.456 

4.482 

2.053 

wmm 

-0.490 

-0.492 

-2.5 

Kim 

5.312 

2.221 

2.228 

-0.589 

-0.591 

-3.0 


6.295 

2.402 

2.409 

-0.697 

-0.700 
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Appendix H 

Formulation for the Problem of a Homogeneous Half Space 
Containing a Crack parallel to the Surface 


The plane strain problem under consideration here is shown in Figure 2.1 with 
/x 2 = /x 1 (7 = 0). Because of the numerical difficulty arising from the formulation 
described in Chapter 2 when the inhomogeneity parameter 7 equals exactly to zero, an 
independent formulation is required for the stability problem of a homogeneous half space 
containing a crack parallel to the surface. 

Substituting 7 = 0 into (2.26)-(2.28), we can write the displacement perturbations 
and v* 2 as 

1 poo 4 

«;(*,») = £C t (a)e"*V“<fa, 

0 < y < h, (HI) 

1 r°° _ - 

v* 2 (x, y) = — Yl mk ^ Ck {oL)z nky z iax da, 

J- oo k=1 


where C 4 are unknown functions of a, rii and m*, i = 1,..., 4, are 

n\ = — A^a], n 2 = - AJ|a|, n 3 = Xl\a\, n 4 = AI|a|, 


m l = i/3 XI, 


i/3 

m2 = A*’ 


m 3 = ~ipX* 2 , 


777-4 = ~ 


ai’ 





(H2) 


and Ai and A 2 are defined by (2.29). The stress perturbations can be obtained by 
substituting (H2) into (2.32) giving 
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<7i„(x, y) = — H J [(1 + a; 2 ) (c 1 e- 1 ;l«l» + C 3 e A SM») 

+ 2(c' 2 e- A i‘l Q ^ + C 4 e A ‘' |Q ^ 


(H3) 


i at' QX da, 



To solve the unknowns functions A\, Ao (in the displacement perturbations for the part 
y < 0) and C\, first we define the density functions (same as in Chapter 2) as 

[®2 (*» + °) “ v i( x > “ 0)] = /i(x), 

— oo < x < oo. (H5a,b) 

— [m;(x, + 0)-<(x, -0)] = / 2 (x), 


Now, by substituting (2.22), (2.28), (2.33) and (H2)-(H4) into (H5) and the following 
boundary and continuity conditions: 

^ 2 y y (Xj h ) = 0, o* 2xy {x, h) = 0, - oo < x < oo, (H6a,b) 

^lyyip^i *^2yy(. X i "F 0), G^ X y{x^ 0) &2xy(> X i "F 0), OO <! X < OO, (H7a,b) 

we can express the six unknown functions A\, A 2 and C 4 in terms of the density 

functions /i and / 2 . From (H5) and the displacement continuity along y = 0, i.e., 

Ui(x, — 0) = u 2 (x, + 0), v\(x, — 0) = i» 2 (x, + 0), |x| > a, (H8a,b) 

we can write the single-valuedness condition as 

fi(x) = 0, |x| > a; f fi(t)dt = 0, * = 1,2. (H9) 

J —a 
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Also, from boundary condition on the crack surface we may write 

0\ yy {x, ” 0 ) = ^ [/ /i (*) dt J Bn ( a ) i /3e la{x ~ t] da 

/ a poo 

f 2 (t)dt B 12 (a)^ x ~^da 

a J -oo 

= 0, — a < x < a, 




, -0)= ^ jT/i (*)*/” Ba(a)e i * < *-‘ ) do 

/ a poo 

f 2 (t)dt I B 22 (a)ipe ia{x - t) dc 

•a J —oo 


0, — a < x < a, 


where 


B] io = 


4A|A£-(1+A| 2 ) 2 
2A$(1-AJ 2 ) ’ 

2\\ [ 4 AIA 2 + (1 + Af ) 2 


-Sl12 = 


B-lO-l — 


(l-Af)[4AiA5-(l + A; 2 ) 2 ]’ 

(1 + Ax 2 ) 2 [4Ai A 2 + (1 + Ai 2 ) 2 ' 
2A^(l-Af)[4AlA^-(l + At 2 )2]’ 

8At(l + Af) 2 _ 

(l-A| 2 )[4AtA* 2 -(l+Af) 2 ]’ 

(1 + AJ 2 ) 4AJA£ + (1 + At 2 ) 2 
(1 - A f) 4\l\* 2 - (1 + A f) 2 


(H10) 


(Hii) 


Bij(a) = BijoSij + B ifl e- 2X ^ h + B ij2 c~ 2 ^ h + B ij3 e~^ +X ^ h , (H12) 

i, j = 1, 2, 


(HI 3) 


(HI 4) 


(HI 5) 


(H16) 


(H17) 


B i22 = Bni, 


(HI 8) 


-Si 23 — — 25 j21, 

(HI 9) 

£211 = — £l21, 

(H20) 

£212 = - £l22, 

(H21) 

£213 = — £l23, 

(H22) 

4A;a;-(i+ a; 2 ) 2 

220 2At(l - AJ 2 ) ’ 

(H23) 

(1 + A! 2 ) 2 [4AjA 2 + (1 + Xf) 2 

2A;(1 - xf) 4X\X* 2 - (1 + Xf) 2 ] ’ 

(H24) 


2A 2 [AXlX* 2 + (l + Xf) 2 \ 

(l-At 2 )[4AlA^ 2 -(l + At 2 ) 2 ]’ 


(H25) 


£22.3 = 


8a;(i + a; 2 ) 2 


(1 — Aj 2 ) [4A|A 2 - (1 + AJ 2 ) 2 ' 


(H26) 


The integral equations (H10) and (Hll) can further be simplified as follows: 

/ q 2 r c 

£ + **(*>*) fib)** = °> 

i — 1,2, — a < x < a. 


where 


u ( T f\ _ Em ^ Et 1 EmE El 

ij{ ’ ; " 4X*fh 2 + {t- xf + 4A fh 2 + (t- x) 2 


(H27) 


+ _ Ml ~*> _ ii = 12 

(A; + A;) 2 fc 2 + ((-*) 2 ’ ’ 3 • 


(H28) 


It can be seen that the singular integral equations (H27) have the same form as (2.66). It 
can further be shown that the coefficients for the Cauchy kernels are the same as in (H27) 
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and (2.26). Also we have k^\ (x, t) = — & 12 OM) which is the same as in the case of 
inhomogeneous coating bonded to homogeneous substrate. To solve the integral 
equations (H9) and (H27) we follow the procedure described in Chapter 3. Following 
Section 3.3.2, a pilot case for the problem of homogeneous half space containing a crack 
parallel to the surface is analyzed to examine the convergence. The computed instability 
strain (eo) CT as a function of approximation terms n is given in Table H.l. Again, it is 
observed that the calculated results by using 8 approximation terms agree very well with 
higher order approximations and therefore, n = 8 is used in this study. Note that because 
the Fredholm kernels in this case have closed-form expressions (i.e., (H28)), the numerical 
computation is greatly simplified. 


Table H. 1 Effect of the number of approximation terms on the instability point 
and phase angle for the problem described in Figure 2.1 with 'yh = 0, is = 0.3, 
and h/a = 0.05. 


Approximation terms n 

Instability load (eo) CT 

Phase angle ip (degree) 

2 

0.001916 

-25.52 

4 

0.001906 

-37.30 

8 

0.001921 

-39.12 

16 

0.001921 

-39.07 

32 

0.001921 

-39.07 
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Appendix I 

The Plate Model Analysis 


1.1 Stability Analysis 



Figure 1.1 Schematic of the plate model for the problem described in Figure 2.1. 


3 






-**-J 



Figure 1.2 Configuration of the column with transformed cross-section. 


To study the buckling instability problem described in Figure 2.1, we may approximate 
the delaminated coating as a plate with built-in ends under fixed-grip compressive strain eo 
as shown in Figure 1.1. The plane strain problem can further be reduce to a FGM column 
of unit width with built-in ends. The elastic modulus of the column can be expressed as 

E'(y) = E[t«, El = (II) 

where E\ is the elastic modulus of the substrate and E\ = 2(1 + v)n\. The problem may 
then be solved by using the Euler column approach. The material inhomogeneity described 
in Equation (II) is considered by using the transformed-section method [88], that is 
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replacing the inhomogeneous column with rectangular cross-section by a homogeneous 
column with stretched width in the cross-section as shown in Figure 1.2. Instead of unit 
width, the transformed cross-section has a width of E*(y)/E{ ( = e 73 '). The position of 
the neutral axis for the transformed section for the case of 7 ^ 0 can be obtained as 


yt^dy 


e^dy 


jhe yh — eJ h + 1 
7(e^ - 1) ’ 


and the moment of inertia for the transformed section is given by 


It = [ (y~ yfdA = [ (y- yft^dy 
Jo Jo 


(e lh - l ) 2 - ( 7 /QV a 

7 3 (e'i' /l - 1 ) 


The axial compressive load for the column corresponding to the fixed-grip compressive 


strain €q can be expressed as 


P= E*e 0 dA 


= /\ 


e™e 0 dy 


E\ e 0 (e 7ft - 1) 


The Euler's buckling load for the column described in Figure 1.2 is given by [87] 


P CT = 


47 X 2 ElIy 
(2a) 2 * 


Substituting (13) and (14) into (15), we have 


E\{€< h - 1 ) 4t x 2 E{ [(e yh - l ) 2 - ( 7 h) 2 & h ] 

7 (6o)cr “ 4a 2 7 3 (e^ - 1 ) 


(c°)ct 7T 2 ( a ) [^ 2 ( ejh+e - lh _ 2 )]> ^ 0 • 


It can be seen that the critical strain (eo) CT corresponding to the buckling instability is an 
even function of the inhomogeneity parameter 7 h. 
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For the case of 7 = 0 which corresponding to a homogeneous coating, one may repeat 
the procedure and obtain the counterpart of (16) as 

At: 2 El /h 3 \ 


o)« - 


4o 2 V12 


(s)- 


(18) 


The critical strain for 7 = 0 can be reduce from (18) and is given by 

7T 2 / h \ 2 




(19) 


Note that (19) may be recovered by taking limits (7 -4 0) of (17). Also it may be seen from 
(17) and (19) that for a fixed value of hf a the maximum of (eo)cr occurs at 7 = 0. 


1.2 Postbuckling Analysis 

By following Hutchinson and Suo [41], the crack problem associated with the 
postbuckling of coating is obtained by coupling the von Karman plate solution and the 
split-beam approach for the fracture mechanics parameters. It is assumed that the ratio of 
crack length (or radius) over coating thickness is very large and the coating thickness is 
much smaller than the one of the substrate. The multi-layered TBC system containing an 
interface crack subjected to thermal or mechanical compressive stress is approximated by 
a coating bonded to a half-space with an interface crack as shown in Figure 1.3(a). The 
loading N and bending moment M on the buckled coating as shown in Figure 1 . 3 (b) is 
then calculated by using the von Karman nonlinear plate formulation. The plane strain 
problem and axisymmetric problem will be discussed separately. 


1.2.1 The Plane Strain Problem 

Again, the inhomogeneity of the coating is considered by using the transformed- 
section method as described in Appendix 1.1. The thermomechanical properties E, u, and 
a are functions of y. The position of the neutral axis and the bending stiffness are given by 
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M 



Figure 1.3 Geometry of the one dimensional blister, (a) the buckled coating, 
(b) local loading of interface crack. 


i 


^ p 



^- 

/ n. y 

1 

1 

7*7 X 


x = -a 


x = a 


Figure 1.4 Configuration of the plate for the postbuckling analysis. 
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yE*{y)dy 
k E* (y)dy 


e*(v) 


1 — v 2 {y)' 


D= / E*(y)(y — y) 2 dy. 
Jo 


The one-dimensional deformation of the plate in Figure 1.4 is characterized by the y- 
displacements, v(x), which is zero in the unbuckled state with pre-stress corresponding to 
the compressive load. From the moment equilibrium of the column, we may write the 


governing equation as 


D 0 + A *S = 0 ’ = 


where AN is the net x-direction resultant compressive force for the buckled column. 
Solving (112) with the following boundary condition: 


we obtain the deflection and the corresponding buckling load (lowest eigenvalue) as 


v(x) = |^( 1+cos ~)’ 


AN = P CT = 


cr — o 

a 1 


respectively, where £ = v(0)/h is the unknown "amplitude". The amplitude £ is 
determined by the condition that the membrane stress in the buckled coating is the same as 
the buckling stress. This condition leads to 

J_ e xx dx = 2a{?^^-), \ = J*E*{y)dy, (115) 


where the nonlinear strain e x is given by 
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Figure 1.5 Geometry of the multi-layered infinite plate. 


_ 1 (dv\ 2 

£ ““ 2 \ii) ' 

By substituting (114) and (116) into (115) we obtain 


\ D i 

(P 41 

lx' 

^ Per VJ 


The bending moment at x = a and the change in resultant force in re-direction can be 
expressed as 


M = D- 


7 i 2 Dh 


t, N = P-P a . 


By substituting (118) into the general split-beam solution given by Hutchinson and Suo 
[41] we can obtain the strain energy release rate as 




>+ 4 )' 


The only undetermined quantity in (117) and (119) is the initial in-plane compressive 
load P. To calculate P, we assumed that initially the crack is closed and there is no shear 
stress in the medium. Then P is approximated by the in-plane compressive load in the 
coating of a multi-layered medium without interface crack as shown in Figure 1.5. It may 
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be seen from Figure 1.5 that for the medium we have — () = —() = 0, and that from 

ox oz 

the compatibility equations, we have 


d%x 

dy 2 


dy 2 


= 0. 


( 120 ) 


Also it may be seen that the out-of-plane stress 

Oyy = 0 . 


( 121 ) 


Solving (120), we obtain 


€ X x — A\y + A2, 


£zz — B\y + B2, 


( 122 ) 


where A\, A 2 , B\, and B 2 are constants and we have 

B\ = B 2 — 0, for plane strain condition. (123) 

The remaining unknown constants A\ and A 2 are determined from the boundary 
conditions associated with the following loading conditions: 

(i) Fixed-grip Loading 

For the case of fixed-grip loading shown in Figure 2.1, we simply have 

Ai — 0, A 2 — — € 0 , for fixed-grip loading. (124) 

By substituting (I21)-(I24) into the generalized Hooke's law, we obtain the stress for the 
coating as 

_ _ E{y) 

l- u 2( y ) e °- (I25) 

The compressive resultant force may be expressed as 

P= ~ f o xx dy = Ae 0 . (126) 

Jo 
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By substituting (126) into (117) and (119) we may then obtain the crack opening and strain 
energy release rate. 

(ii) Temperature Loading 

By substituting (I21)-(I23) into the generalized Hooke's law containing temperature 
effects, we can obtain the thermal residual stress for the layered system shown in Figure 
1.5 as 

° ixx = 1 -l?(y) ^ AlV + ~ y { . x <y< y u (127) 

where the subscript i denotes the i-th layer and 

n 

2/o = 0, yi = Y^hi. (128) 

i=l 


From the force and moment equilibrium on each cross-section, i.e., 



t =i Jyi -1 


tf: 


(tixxydy = 0 , 


(129) 


we can calculate the unknown constants A\ and Ao and, in turn, determine the stresses 
Oi xx . Assuming that there are m layers above the interface crack, the resultant 
compressive load can be expressed as 

_!L rvi 

p= Y, ( I3 °) 

i—n—m Jyi-\ 

The results for the postbuckling regime can be then obtained by substituting (130) into 
(117) and (119). 
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1.2.2 The Axisymmetric Problem 


z 


A 



Figure 1.6 Configuration of the circular plate for the postbuckling analysis. 


The problem of an axisymmetric blister is investigated by following the same 
procedure used in the plane strain case. In this study only the interface crack problem for a 
homogeneous coating is discussed. From the moment and force equilibrium of the circular 
plate described in Figure 1.6 we have 


ri (d 3 w i 1 d 2 w 

1 dw^ 

, dw 



V dr 3 * r dr 2 

r 2 dr > 

i + ajv '* =0, 

AN t = P-N t , 

( 131 ) 



Ehr / dw\ 2 



— 1 
dr 

( r V) 

II 

3 

<N 

+ 

0 . 

( 132 ) 


The boundary and symmetric conditions are given as follows: 

dw n , . dw 

*„o = °’ "’ (a) = 0, * =°’ 


(I33a) 


r r d i dN 

=0, = 0, (I33b) 


For nontrivial solution for Equation (131) under the boundary conditions (I33a), we have 

(A^ir)cr — Per ~ a 2fr2 ’ (134) 

where (i = 3.8317 is the first nontrivial zero of Ji (x), which is the Bessel function of the 
first kind of order one. The associated buckling mode is 
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(135) 


wi (r) 


J 0 (Atr/g) - Jo(jtt) - 

1 - J 0 (At) 


h. 


where J 0 is the Bessel function of the first kind of order zero. Note that the buckling mode 
is normalized such that W\ (0) = h. To examine the initial postbuckling behavior a 
perturbation analysis is performed. It is assumed that 

w(r) = (r) + fw 2 (r) + fw 3 (r) + 0(£ 4 ), (136) 

N t = £N n + £ 2 )V r2 + 0(£ 3 ). (137) 

Substituting (136) and (137) into (131) and (132) and writing them as power series in £ we 
can equate the coefficients to zero to generate an ordered series of equilibrium equations 
and boundary conditions. Solving these equations we obtain [41] 

N a = 0, N i2 = CjPcr, (138) 


where c\ = 0.2473(1 + v) -\- 0.2231(1 — u) 2 . From (137) and (138) we may write 

. 6 r 1 / a N] 1 / 2 

f-s-bb- 1 )] • (I39) 


Again, substituting the following quantities into split-beam solution [41] 

. d 2 w I 


M = D 


dr 2 


, N = N T \ r=a , 


(140) 


we obtain the asymptotic strain energy release rate for initial postbuckling as 



1 

1 + 0.9021(1 — v) ’ 


(141) 


where Gq = (1 - u)P 2 /Eh. 

In the axisymmetric problem, only the compressive stress induced by temperature 
change will be discussed here. To determine the initial in-plane compressive load, we 
follow the procedure described in previous section. First we calculate the stresses in 
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Cartesian coordinate and transform them into cylindrical coordinate to obtain <j rr . For the 
medium described in Figure 1.5 under "axisymmetric" condition, we may write 

®yy ~ Oj txx = ^zz ~ A\y + A 2 , (142) 

where A\ and A 2 are the unknown constants. Substituting (142) into the generalized 
Hooke's law, we obtain 

<?ixx = cr izz = El ^ (Aiy + A 2 ) - - E ( r ai(y)AT, y { - x <y<y it (143) 
i -Vi(y) i -My) 

where y* is defined in (128). The unknowns A\ and A 2 can be obtained by substituting 
(143) into (129) and solving the system of equations. Under the stress state described by 
(142) and (143) we can simply write 

rr = ®ixx > ® = lj (144) 

The resultant compressive load for the homogeneous coating above the crack (layer n) 
can then be simply expressed as 

P = cr nxx h. (145) 

The results for the postbuckling regime can then be obtained by substituting (145) into 
(139) and (141). 
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